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previous  theory  are  briefly  discussed.  Typical  driver  results 
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representation  of  any  nuclear  weapon  neutron  energy  spectrum 
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1.  INTRODUCTION 


Various  Army  electronic  systems  may  be  exposed  to  nuclear  explosion 
effects  in  a  tactical  battlefield.  Included  among  these  nuclear  effects 
is  the  nuclear  electromagnetic  pulse  (EMP),  which  is  a  transient  broad¬ 
band  electromagnetic  field  capable  of  damaging  or  upsetting  electronic 
equipment.  To  predict  the  signature  of  the  BMP  generated  by  a  nuclear 
burst,  it  is  necessary  to  determine  the  physical  parameters  that  induce 
the  EMP.  These  parameters,  called  EMP  drivers,  are  time  and  space 
varying  ionization  and  Compton  electron  currents  in  the  nuclear  radia¬ 
tion  field  around  the  burst.  Once  the  ionization  and  the  currents  are 
specified,  it  is  usually  possible  to  solve  some  form  of  Maxwell's  equa¬ 
tions  for  the  EMP  generated.  This  report  summarizes  the  development  of 
a  software  package  specifying  EMP  drivers  produced  by  neutrons  emitted 
from  a  near-surface  burst  in  the  air.  The  package  was  designed  for  use 
with  the  NEMP  computer  code  used  at  the  Harry  Diamond  Laboratories 
(HDL).  A  separate  report1  summarizes  the  development  of  a  similar 
software  package  for  EMP  drivers  due  to  prompt  gamma  radiation. 


2.  MOTIVATION 

The  DIP  drivers  are  due  to  one  component  arising  from  prompt  gamma 
radiation  and  a  second  component  arising  from  neutron  radiation.  The 
neutron  component  includes  elastic  recoil  ionization  (frequently  called 
"heating")  by  air  constituent  nuclei  and  also  includes  effects  of 
secondary  gamma  radiation  produced  by  a  number  of  neutron  capture  and 
inelastic  collision  reactions  in  the  air  and  the  ground.  The  develop¬ 
ment  of  H4P  drivers  by  secondary  gamma  rays  is  physically  similar  to  the 
development  for  prompt  gamma  radiation,  although  the  magnitude,  the 
direction,  and  the  time  dependence  are  markedly  different. 

Maxwell's  equations  for  EMP  from  a  near-surface  burst  are  solved  by 
the  NEMP  computer  code . 2 ' 3  The  NEMP  code  requires  EMP  drivers  to  be 
specified  in  a  volume  extending  to  several  kilometers  from  the  burst 
point  and  for  times  extending  to  many  milliseconds  after  the  instant  of 
burst.  The  EMP  drivers  must  be  incorporated  into  the  code  as  smooth 
fits  to  results  of  quasi-analytic  or  Monte  Carlo  predictions.  The  NEMP 

1 William  T.  Wyatt,  Jr.,  A  Near-Surface  Burst  EMP  Driver  Package  for 
Prompt  Gamma-Induced  Sources,  Harry  Diamond  Laboratories  HDL-TR-1931 
(1980), 

2H.  J.  Dongle y,  C.  L.  I ongmire,  and  K .  S.  Smith,  Development  of  NEMP 
(U),  Mission  Research  Oorp,,  Santa  Barbara,  C A,  HDL-CR-7 5-00 1-1  (April 
1975).  (SECRET— RESTRICTED  DATA) 

3H.  J.  long  ley  and  K.  S.  Smith,  Developments  in  NEMP  for  1977  ( U ), 
Mission  Research  Oorp.,  Santa  Barbara,  CA,  HDL-CR-77-0022-1  (January 
1978).  ( SECRET— RESTRICTED  DATA) 


code  in  many  respects  is  an  extension  of  the  LEMP  computer  code**  for  EMP 
from  surface  bursts,  which  was  initially  developed  from  about  1966  to 
about  1969.  The  BMP  drivers  in  the  earlier  LEMP  code  could  have  been 
modified  and  used  in  the  NEMP  code  with  certain  limitations.  These 
limitations  are  discussed  briefly,  and  reasons  are  given  for  the  impor¬ 
tance  of  the  new  results  reported  here. 

In  this  study,  the  neutron- induced  EMP  drivers  are  divided  into  five 
categories:  (1)  neutron  elastic  (recoil)  ionization,  (2)  fast  neutron 

air  inelastic  collision  (n,n"Y)  drivers  (both  Compton  electron  currents 
and  ionization),  (3)  fast  neutron  ground  inelastic  collision  (n,n'Y)  and 
capture  (n,y)  drivers,  (4)  thermal  neutron  ground  capture  (n,Y)  drivers, 
and  (5)  thermal  neutron  air  capture  (n,y)  drivers.  Of  these,  one  would 
expect  the  NEMP  geometry  to  introduce  complications  to  LEMP  code 
prescriptions  for  categories  (1),  (2),  and  (4),  since  the  proximity  of 
the  ground  is  different.  Category  (5)  was  not  treated  originally  by  the 
LEMP  code.  Category  (3)  would  be  quite  different  for  the  NEMP  code  and 
probably  of  reduced  magnitude.  However,  the  cross  sections  for  most  of 
the  above  reactions  have  been  subjected  to  concentrated  study  and  revi¬ 
sion  with  the  advent  of  the  Evaluated  Nuclear  Data  Pile  (ENDF-B, 
distributed  Brookhaven  National  Laboratory,  Brookhaven,  NY).  It  is 
likely  that  new  neutron  transport  studies  based  on  recent  ENDF-B  cross 
sections  would  produce  results  different  from  those  upon  which  the  LEMP 
code  neutron- induced  HIP  drivers  were  based.  Further,  by  omission  of 
the  thermal  neutron  air  capture  (n,Y)  drivers,  the  LEMP  drivers  were 
limited  to  times  less  than  about  1  ms  after  the  burst. 

The  LEMP  authors4 5  developed  a  newer  and  much  more  complete  prescrip¬ 
tion  of  categories  (1)  through  (5)  based  on  results  of  Sargis  and 
others.5  Factors  weighing  against  the  use  of  these  prescriptions  in  the 
NEMP  code  include  the  following:  First,  a  newer  revision  (Bound  3)  of 
ENDF-B  was  released  in  1973.  Second,  the  only  burst  heights  considered 
by  Sargis  were  0,  200,  and  500  m,  whereas  the  NEMP  code  operates  princi¬ 
pally  at  tactical  burst  heights  between  0  and  200  m.  Third,  the  new 
LEMP  EMP  driver  curve  fits  were  considered  by  this  writer  to  neglect 
certain  important  ground-air  interface  effects  such  as  time-dependent 
depletion  and  enhancement.  Fourth,  the  transport  results5  were  based  on 

4H.  J.  bongley  and  C.  L.  bongmire,  Development  and  Testing  of  LEMP  1, 
Los  Alamos  Scientific  Laboratory ,  MM,  LA-4346  (April  1970)% 

SH.  J.  bongley,  C,  L.  Longmire,  J,  S.  Malik,  R,  M.  Hamilton,  R.  N. 
Marks,  and  K.  S,  Smith,  Development  and  Testing  of  LEMP  2,  A  Surface 
Burst  EMP  Code  (U),  Mission  Research  Oorp,,  Santa  Barbara,  CA,  DMA  4097T 
(December  1976),  (CONFIDENTIAL) 

5D.  A,  Sargis,  E,  R,  Parkinson,  J,  N.  Mood,  R,  E,  Dietz,  and  C ,  A, 
Stevens,  Late-Time  Sources  for  Close-In  EMP,  Science  Applications,  Inc,, 
La  Jolla,  CA,  DMA  3064F  (August  1972), 


"energy  grouped"  cross  sections#  whereas  a  "point"  cross-section  treat¬ 
ment  should  be  more  accurate  at  deep  penetrations*  Fifth,  a  single 
typical  source  neutron  energy  spectrum  was  used  in  the  I£MP  drivers, 
with  a  parameter  available  for  adjusting  the  concentration  of  high- 
energy  (14-MeV)  neutrons;  a  more  desirable  package  would  allow  flexible 
representation  of  arbitrary  source  neutron  energy  spectra* 

For  these  reasons,  new  Monte  Carlo  transport  calculations  have  been 
done  to  describe  more  accurately  the  neutron- induced  Q(P  drivers. ^  The 
extent  of  these  calculations  is  summarized  in  table  1.  Edited  results 
of  these  calculations  were  transmitted  to  this  writer  on  computer  magne¬ 
tic  tape  and  used  to  derive  the  EMP  driver  software  package  reported 
here. 

TABLE  1.  MONTE  CARLO  DATA  CHARACTERISTICS 


Source  neutron 

Tine  bin 

Burst 

EHP 

driver 

energy  band 

(Mav) 

boundary 

<«> 

height 

(a) 

0.0  to  0. It 

0.0  to  1.01-7)* 

1 

Neutron  direct  ionisation 

O.tt  to 

0.55 

1.0(-7)  to  2.151-7) 

50 

High-energy  air,  radial  current 

0*55  to 

1.11 

2.151-7)  to  4.641-7) 

100 

High-energy  air.  theta  current 

1.11  to 

1.83 

4.641-7)  to  1.01-6) 

200 

High-energy  air,  ionization 

1.83  to 

2.35 

1.01-6)  to  2.151-6) 

Low-energy  air,  radial  current 

2.35  to 

4.07 

2.151-6)  to  4.641-6) 

Low-energy  air,  theta  current 

4.07  to 

6.36 

4.641-6)  to  1.01-5) 

Low-energy  air,  ionization 

6.36  to 

8.19 

1.01-5)  to  2.151-5) 

High-energy  ground,  radial  current 

8.19  to 

15.0 

2.151-5)  to  4.641-5) 

High-energy  ground,  theta  current 

4.641-5)  to  1.01-4) 

High-energy  ground,  ionization 

1.01-4)  to  2.151-4) 

Low-energy  ground,  radial  current 

2.151-4)  to  4.641-4) 

Low-energy  ground,  theta  current 

4.641-4)  to  1.01-3) 

Low-energy  ground,  ionization 

1.01-3)  to  2.151-3) 

2.151-3)  to  4.641-3) 

4.641-3)  to  1.01-2) 

1.01-2)  to  2.151-2) 

2.151-2)  to  4.641-2) 

4.641-2)  to  1.01-1) 

*Rma i  am  1.0  x  10~7. 

'HT.  S,  Schechter  and  M.  0.  Cohen,  Energy  Deposition  Rates  and  Compton 
Electron  Currents  from  Low-Altitude  Bursts  as  a  Function  of  Source 
Energy,  Mathematical  Applications  Group,  Inc.,  El  ms ford,  NY,  HDL-CR-77 - 
020-1  (November  1977 )• 
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3.  APPROACH 


The  edited  results  of  Monte  Carlo  calculations  contain  statistical 
fluctuations  inherent  in  the  Monte  Carlo  method.  The  NEMP  code  demands 
smooth  prescriptions  of  the  EMP  drivers  to  obtain  useful  results.  One 
may  attempt  to  fit  smooth  general  functions  (such  as  polynomials)  to  the 
Monte  Carlo  results  by  a  least-squares  fitting  technique,  for  example, 
or  one  may  attempt  to  fit  the  Monte  Carlo  results  with  certain  appro¬ 
priate  smooth  functional  forms.  By  using  a  polynomial  of  sufficiently 
high  degree,  the  former  technique  may  allow  the  Monte  Carlo  results  to 
be  fitted  more  faithfully  than  the  latter,  but  the  closer  fit  may  mean 
that  the  "noise"  is  being  fitted  instead  of  the  underlying  smooth  "true 

answer."  The  same  difficulty  can  arise  by  using  the  latter  technique  if 

too  many  degrees  of  freedom  are  permitted.  On  the  other  hand,  if  the 

polynomial  degree  is  not  high  enough,  the  fit  may  be  a  poor 

representation  of  the  true  answer.  In  general,  however,  using 

appropriate  smooth  functional  forms  will  involve  far  fewer  degrees  of 
freedom  than  using  general  functions  if  the  functional  forms  jure  well 
chosen  and  if  the  fit  is  a  global  fit  over  the  entire  space-time  volume 
of  interest  rather  them  a  piecewise  fit  over  many  smaller  subdivisions 
of  the  space-time  volume. 

The  strength  of  the  technique  of  using  well-chosen  functional  forms 
is  that  physical  variables  usually  do  behave  according  to  relatively 
simple  physical  principles  that  are  approximately  described  by  rela¬ 
tively  simple  functional  forms.  A  global  fit  of  properly  chosen  func¬ 
tional  forms  to  Monte  Carlo  results  could  reduce  the  error  below  the 
statistical  error  of  the  Monte  Carlo  calculation.  The  fitting  algorithm 
should  heavily  might  low-variance  results  over  high-variance  results  to 
reduce  the  error. 

The  new  edited  Monte  Carlo  results  for  neutron- induced  EMP  drivers 
consist  of  three  physical  quantities  (energy  deposition  and  radial  and 
polar  components  of  Compton  electron  current)  defined  for  19  time  inter¬ 
vals  within  63  spatial  volume  detectors,  for  five  categories  of  reaction 
(recoil,  air  inelastic,  etc.),  four  burst  heights,  and  nine  source 
neutron  energy  bands  (table  1).  There  are  646,380  quantities  to  be 
fitted.  Apparently,  the  fitting  algorithm  must  be  highly  automated  to 
allow  useful  results  to  be  obtained  with  a  realistic  investment  of  human 
effort. 


4.  CURVE-PITTING  PROCEDURE 

Thirteen  distinct  functions  describe  the  time  and  space  dependence 
of  13  EMP  drivers: 


( 1 )  Ionization  in  the  air  from  neutron  elastic  collisions  and 
charged  particle  production  (such  as  (n,p)  and  (n,a)  reactions  that 
deposit  kinetic  energy  locally) 

(2)  Radial  Compton  electron  current  in  the  air  (measured  in  the 

polar  coordinate  system  centered  at  the  burst)  due  to  high-energy  neu¬ 
tron  reactions  in  the  air  (neutron  energy  more  than  0. 1  MeV) 

(3)  Polar  Compton  electron  current  in  the  air  for  the  same  as  (2) 

(4)  Ionization  in  the  air  for  the  same  as  (2) 

(5)  Radial  Compton  electron  current  in  the  air  due  to  low-energy 
neutron  reactions  in  the  air  (neutron  energy  less  than  0.1  MeV) 

(6)  Polar  Compton  electron  current  in  the  air  for  the  same  as  (5) 

(7)  Ionization  in  the  air  for  the  same  as  (5) 

(8)  Radial  Compton  electron  current  in  the  air  (measured  in  the 

polar  coordinate  system  centered  on  the  ground  beneath  the  burst)  due  to 
high-energy  neutron  reactions  in  the  ground 

(9)  Polar  Compton  electron  current  in  the  air  for  the  same  as  (8) 

(10)  Ionization  in  the  air  for  the  same  as  (8) 

(11)  Radial  Compton  electron  current  in  the  air  due  to  low-energy 
neutron  reactions  in  the  ground 

(12)  Polar  Compton  electron  current  in  the  air  for  the  same  as  (11) 

(13)  Ionization  in  the  air  for  the  same  as  (11) 

Drivers  (2),  (3),  and  (4)  are  due  principally  to  high-energy  inelas¬ 
tic  collisions  with  atmospheric  nitrogen.  Drivers  (5),  (6),  and  (7)  are 
due  to  thermal  capture  by  atmospheric  nitrogen.  Drivers  (8),  (9),  and 
(10)  are  due  to  inelastic  collisions  with  the  ground  and  fast  neutron 
capture  by  the  ground.  Drivers  (11),  (12),  and  (13)  are  due  to  thermal 
neutron  capture  by  the  ground. 

A  mathematical  notation  is  used  to  simplify  the  succeeding  discus¬ 
sion.  The  vector  functional  F  represents  the  drivers.  F  is  dependent 
on  time  (t)  and  space  (x,y,z).  Subscript  d  running  from  1  to  13  denotes 
which  of  the  drivers  is  represented: 

Fd(t,x,y,z) 
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Time  t  is  the  local  time,  retarded  by  the  speed  of  light  from  the 
coordinate  center. 


The  functional  Fd  is  a  functional  form  incorporating  a  certain 
number,  nd,  of  parameters  p.  This  is  written  as 

Fd(t,x,y,z,Pdi)  ' 

where  pdi ,  i  =  1  to  nd,  is  interpreted  as  the  parameter  vector  of  length 
nd  for  driver  d.  In  general,  nd  is  not  the  same  for  different  d,  but  is 
approximately  10  in  this  report. 

Parameters  pdi  are  evidently  dependent  on  the  energy  of  the  source 
neutrons  producing  the  drivers  and  on  the  height  of  burst.  Other  varia¬ 
bles  are  ignored,  such  as  variations  in  water  content  of  the  air  and  the 
ground  and  variations  in  ground  composition.  The  parameters  may  then  be 
written 

Pdi(j*h) 

or 

Pdij(h)  * 

where  j  denotes  the  ordinal  of  the  source  neutron  discrete  energy  band 
and  h  is  the  height  of  burst.  In  this  way,  p  is  a  discrete  function  of 
the  source  neutron  energy  (band)  and  a  continuous  function  of  the  height 
of  burst.  It  thus  appears  that  the  EMP  driver  fits 

Fd(t,x,y,zfpdij) 

are  defined  when  parameters  Pd^j  are  evaluated  for  a  specific  height  of 
burst, 

pdij  ”  pdij^  * 

The  total  number  of  functions  Pdij(h)  to  be  determined  is 
13  nd  9  13  nd 

£  E  £  i  =  £  £  ®  • 

d«1  i=1  j-1  d= 1  i=1 


Since  nd  ~  10,  this  total  is  roughly  1100,  Thus,  there  are  over  a 
thousand  such  functions  to  be  specified.  Determination  of  this  many 
functions  was  a  formidable  undertaking  and  required  considerable 
computer  time  (about  40  hr  on  an  IBM  System/370  Model  168), 
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The  following  strategy  was  adopted: 


a.  Through  previous  experience  supported  by  trial  and  error,  13 
functional  forms  were  constructed  to  constitute  driver  functionals  Fd. 
These  were  functions  of  time  and  space  and  were  based  on  a  certain 
number,  nd,  of  parameters. 

b.  Monte  Carlo  predictions  for  the  13  drivers  were  available  as 
functions  of  time  for  each  of  63  spatial  volune  detectors  around  the 
burst  for  each  combination  of  four  burst  heights  and  nine  source  neutron 
energy  bands.  For  a  particular  source  neutron  energy  band,  each  driver 
functional  was  fitted  to  the  time  and  space  dependence  of  the  Monte 
Carlo  data,  for  each  of  the  four  sets  of  data  corresponding  to  the  four 
heights  of  burst  (1,  50,  100,  and  200  m) .  The  fit  was  obtained  through 
an  optimization  process  by  adjustment  of  each  of  the  nd  parameters  for 
the  driver. 

c.  Four  values  ,  Pdij^hk^»  k  =  1  to  4,  were  thus  derived  for 
each  parameter  Pd^j  corresponding  to  the  four  heights  of  burst.  These 
four  values  were  approximated  by  simple  functions  of  the  height  of 
burst.  In  some  cases,  a  constant  function  was  used;  in  other  cases,  a 
rational  fraction  was  used.  Certain  constraints  were  imposed  on  the 
rational  fraction  to  insure  satisfactory  behavior  of  the  approximation. 

d.  Since  parameter  functions  pdij(h)  only  approximated  values 

Pdij(hk),  It  *  1  to  4,  it  was  desirable"*  to  readjust  the  amplitudes  of 
the  driver  functionals  Fd( t,x,y,z;pdjj (h) )  to  obtain  a  closer  fit  to  the 
Monte  Carlo  data.  Since  the  amplitude  was  by  choice  one  of  the  param¬ 
eters— specifically,  pd1  . ,  this  readjustment  involved  simply  reoptimi¬ 
zing  the  fit  to  the  Monte  Carlo  data  by  adjustment  of  this  one  param¬ 
eter,  Once  this  adjustment  was  done  for  each  of  the  four  heights 

of  burst,  amplitude  parameter  function  pd1 j(h)  was  refitted  to  these 
four  improved  values  for  Pdlj(hk)* 

e.  This  completed  the  definition  of  all  parameter  functions 
Pdia(h)  for  a  particular  neutron  energy  band,  j.  Process  b  to  d  was 
repeated  for  each  of  the  nine  energy  bands. 

This  strategy  was  implemented  through  the  use  of  two  specially 
developed  computer  codes,  BIGFIT  and  HOBFIT.  BIGFIT  performs  a  global 
minimization  (that  is,  simultaneously  over  all  volume  detectors  and  time 
bins)  of  an  error  function  estimating  the  "badness  of  fit"  of  driver 
functionals  Fd  to  the  Monte  Carlo  data.^  HOBFIT  approximates  the  burst 
height  dependence  of  parameter  values  Pdijth^)  with  either  a  constant 
function  or  a  rational  function.  BIGFIT  performs  the  minimization  for 
specified  parameters  of  the  functional.  In  practice,  the  specified 
parameters  were  usually  either  the  amplitude  parameter  only  or  all 
parameters  of  the  functional • 
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Bie  curve-fitting  procedure  is  illustrated  in  figure  1.  This  proce¬ 
dure  was  performed  for  each  source  neutron  energy  band  and  each  driver 
functional,  or  9  x  13s  117  times,  except  where  the  driver  was  exactly 
zero  for  an  energy  band.  (For  example,  no  high-energy  neutron  reactions 
can  occur  for  a  low-energy  source  neutron.)  First,  starting  with  trial 
guesses  for  the  parameter  values,  BIGFIT  was  run  four  times  (once  for 
each  burst  height)  to  fit  the  driver 
functional  to  the  Monte  Carlo  data 
for  that  case.  All  parameters  were 
adjusted  to  obtain  the  fit.  The 
functional  parameter  values  so  de¬ 
termined  were  stored  in  temporary 
data  file  FITFILE.  Second,  HOBFIT 
obtained  functions  approximating  the 
burst  height  dependence  of  the 
parameters  stored  on  FITFILE,  saving 
these  functions  on  permanent  data 
file  HOBFILE  A.  Third,  the  param¬ 
eter  functions  were  evaluated  at  the 
four  burst  heights,  and  the  result¬ 
ing  parameter  values  were  used  as 
BIGFIT  was  rerun  four  times  (once 
for  each  burst  height)  to  fit  the 
driver  functional  to  the  Monte  Carlo 
data,  this  time  adjusting  only  the 
amplitude  parameter  of  the  func¬ 
tional.  The  four  amplitude  param¬ 
eter  values  so  obtained  were  stored  Figure  1#  Fitting  algorithm.  Final 
in  temporary  data  file  AMPFIT.  amplitude  parameter  fits 

Fourth,  HOBFIT  obtained  a  function  are  created  on  nounm 

approximating  the  burst  height  de-  B.  other  parameter  fits 

pendence  of  the  revised  amplitude  are  created  on  hqbfile 

parameter  values  on  AMPFIT,  saving  A.  T^g  algorithn  is 

this  function  on  permanent  data  file  applied  to  each  driver 

HOBFILE  B.  Thus,  the  final  ampli-  for  every  80urce  neutron 

tude  parameter  fits  are  stored  on  energy  band. 

HOBFILE  B,  and  the  other  parameter 
fits  are  stored  on  HOBFILE  A. 


5.  MINIMIZATION  OF  FIT  ERROR 

The  goal  of  the  curve-fitting  procedure  in  this  work  is  to  obtain  a 
global  fit  to  the  Monte  Carlo  data  for  the  EMP  drivers — that  is,  to  fit 
a  single  functional  to  the  data  over  the  entire  space-time  grid  of 
volume  detectors  and  time  bins.  The  global  fit  is  desired  rather  than  a 


piecewise  fit  to  the  space  or  time  dependence  of  the  data.  For  a  global 
fit,  a  penalty  function  is  constructed  that  measures  the  suitably 
weighted  error  or  the  disagreement  between  the  data  and  the  curve  fit. 
The  curve  fit  is  improved  by  reducing  the  magnitude  of  the  penalty 
function  by  adjustment  of  the  fit  parameters.  A  best  fit  is  obtained 
when  the  penalty  function  arrives  at  a  (minimum)  extremun.  The  penalty 
function  is  a  function  of  the  several  parameters  for  the  driver  func¬ 
tional  and,  with  suitable  smoothness,  is  amenable  to  minimization  by 
standard  optimization  algorithms. 

To  illustrate  the  curve-fitting  procedure,  a  specific  example  of  a 
driver  functional  is  examined.  Neglecting  air  density  and  mass-scaling 
effects  for  simplicity  of  presentation,  the  functional  for  ionization  by 
high-energy  neutron  inelastic  collisions  in  the  atmosphere  (that  is, 
"air  inelastic"  gamma  ray  production)  is 

I(x,y,z,t)  =  pjACD 

where 

x,y,z  =  Cartesian  coordinates  of  the  observer  measured  from  the 
coordinate  center  on  the  ground  beneath  the  burst, 

t  *  retarded  time  at  the  point  (x,y,z"), 

Z**  =  z  -  h, 

h  =  burst  height, 

Pj  =  amplitude  parameter, 

A  *  1  -  P2  exp(-zp3), 

C  *  exp(-tp4)/(tPs  +  p6), 

D  *  exp(-rp?)/(r2  +  pg), 

r  *  [x2  +  y2  +  (z')2]  1 ^2. 

This  particular  driver  functional  is  relatively  simple  and  incorpo¬ 
rates  range  attenuation  factor  D,  time  dependence  factor  C,  and  ground 
depletion  factor  A.  These  factors  are  typical  of  other  driver  func¬ 
tionals,  as  well,  although  other  functionals  may  also  include  time- 
dependent  or  range- dependent  ground  depletion  (or  enhancement) 
factors.  Different  time-dependence  factors  are  used  also. 


In  this  functional ,  eight  parameters  are  used  so  that  the  penalty 
function  is  a  function  of  eight  variables*  The  penalty  function  is 
defined  as  function  P  of  parameter  vector  p^,  k  =  1  to  8: 

■<*>■&  ([^-'(vvvvp,)]2^!  • 

where  the  summation  is  performed  over  all  i  and  j  and  where 

i  =  a  volume  detector  index  running  from  1  to  63, 

j  =  a  time  bin  index  running  from  1  to  19, 

fij  =  the  driver  Monte  Carlo  score  in  the  ith  detec¬ 

tor  and  jth  time  bin, 

f  (xi/yi'zi'tj =  the  inte9ral*  over  the  ith  detector  and  jth 
'  time  bin,  of  driver  functional 

F<xi»yi»2i»t j ipk)  (<*  subscript  omitted  from  P) 
evaluated  for  parameter  vector  p^» 

wii  =  a  nonnegative  weight  facter  defined  as 

A 

fij  *  the  total  ionization  score  in  the  ith  detector 

and  the  jth  time  bin  due  to  all  drivers  (for 
ionization,  these  are  high  and  low  energy  n-y 
reactions  in  air  and  in  ground  and  neutron  di¬ 
rect  ionization) , 

wij  *  a  user-defined  nonnegative  weight  function  «n- 

phasizing  the  importance  of  certain  spatial 
regions  (typically,  about  10.0  near  the  ground 
and  decreasing  to  about  0.25  to  1.0  at  1  km 
above  the  ground) , 

vi-4  =  the  variance  of  score  f  4  4 . 


Score  f^j  would  be  reproduced  exactly  by  integral  f(x^,y^,*i,ta»pk) 
if  driver  functional  F  were  an  exact  fit  to  the  driver  and  if  the  vari¬ 
ance  associated  with  Monte  Carlo  score  fjj  were  zero.  In  this  instance, 
P(P]()  would  be  zero.  In  actuality,  functional  F  is  not  an  exact  fit, 
and  variance  V  is  proportional  to  the  square  of  a  fractional  deviation 
averaging  about  20  percent  for  the  data  given. 


It  is  apparent  that  minimization  of  P(pk)  yields  a  weighted  least- 
squares  fit.  The  use  of  the  inverse  of  variance  increases  the 

weight  of  low-variance  Monte  Carlo  data,  which  is  expected  to  be  more 
accurate  than  high-variance  data  in  inverse  proportion  to  the 
variance.  The  use  of  ratio  f? j/ffj  causes  the  relative  importance  of 
this  driver  to  be  considered  in  comparison  with  all  the  other 
contributing  drivers.  If  the  ratio  is  small,  the  weight  is  diminished 
accordingly;  if  the  ratio  is  near  unity  (that  is,  no  other  driver  is 
contributing  significantly) ,  the  weight  is  maximized.  User-defined 
factor  w^  is  used  to  force  the  best  fit  to  be  obtained  near  the  ground 
(large  w^)  and  a  moderately  good  fit  far  from  the  ground  (small  w^) .  A 
very  good  fit  near  the  ground  is  desired  because  most  BfP  targets  (that 
is,  electronic  systems)  are  found  on  or  near  the  ground. 

To  evaluate  the  quadruple  integral. 


f  (xi'yi'*i'tj,pk)  =  ///  <*x  * y  **  J  at  F(xi#yirxi»tj;pk)  , 

V .  1 1  —  1 

i  1 


over  detector  volume  and  time  bin  < fcj-i j  )  *  a  numerical  integration 
must  be  done.  This  numerical  integration  must  be  done  efficiently 
because  integral  f  must  be  evaluated  63  *  19  «  1197  times  for  each 
evaluation  of  penalty  function  P.  Penalty  function  P  must  typically  be 
evaluated  several  hundred  times  to  locate  a  useful  minimum  and  the 
associated  best  fit.  Thus,  the  numerical  integration  may  be  required 
about  10s  times  to  complete  one  of  the  necessary  4  x  13  x  9  =  468  runs 
of  BIGFIT.  This  is  about  107  integral  evaluations  for  the  entire  fit¬ 
ting  problem.  Evaluation  of  functional  F  requires  roughly  10~3  s  of 
computer  central  processor  (CP)  time,  so  that  about  lO^N  seconds  of  CP 
time  are  forecast  if  N  evaluations  of  functional  F  are  needed  to  obtain 
its  numerical  integral  f.  Clearly,  N  must  be  as  near  unity  as  possible 
to  avoid  prohibitive  computer  time  requirements.  Accordingly,  a  three- 
point  Gauss -Legendre  quadrature  was  used  for  the  time  integration,  and  a 
special,  highly  accurate  two-point  Gaussian  quadrature  (with  adjustable 
weight  function)  was  developed  by  this  writer  and  used  for  the  volime 
integration.  Appendix  A  describes  the  two-point  scheme.  Several  other 
techniques  were  used  to  reduce  the  overall  CP  time  requirement  to  about 
40  hr  for  the  entire  curve-fitting  project. 

With  the  penalty  function  thus  defined  for  the  driver  functional,  it 
is  necessary  to  minimize  the  penalty  function  by  adjustment  of  the 
parameter  vector.  A  highly  efficient  minimization  algorithm  is  desired, 
in  the  sense  of  using  the  fewest  possible  evaluations  of  the  penalty 
function.  A  number  of  different  minimization  algorithms  were  tested. 
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The  most  efficient  algorithm  for  this  application  was  determined  to  be  a 
modified  conjugate  gradient  method.  This  writer  adapted  the  well-known 
conjugate  gradient  method®  to  use  finite  differences  to  estimate  the 
gradient  matrix  and  to  cycle  indefinitely  during  a  single  descent  proce¬ 
dure  instead  of  repeating  the  descent  procedure  after  cycling  once  for 
each  element  in  the  vector  of  independent  variables.  A  very  efficient 
algorithm  with  quadratic  convergence  also  was  developed  to  search  for  a 
minimum  along  a  chosen  direction,  in  connection  with  the  descent  proced¬ 
ure.  Appendix  B  contains  a  listing  of  the  FORTRAN  version  of  these 
descent  algorithms.  It  was  found  that  minimization  of  the  penalty 
function  proceeded  much  more  rapidly  if  amplitude  parameter  pi  were 
optimized  analytically  during  evaluation  of  penalty  function  PCp^). 
Optimum  amplitude  a  (to  minimize  P(p^)  with  respect  to  the  amplitude 
parameter)  can  easily  be  shown  to  be  trial  amplitude  afc  multiplied  by 
the  factor 


where,  as  previously  explained, 

f 

f^j  “  the  data  value  for  volume  detector  i  and  time  bin  j, 

f(  )  =  the  corresponding  value  of  the  fit  using  a  nonzero  trial 
value  (completely  arbitrary)  for  p  , 

g^j  »  corresponding  complete  weight  function  W^j/V^j  used  in 
penalty  function  P(pJc) • 

When  p  is  chosen  in  this  way  to  be  a  =  K and  when  the  minimization  of 
P(pk)  is  done  for  the  other  pk  exclusive  of  p^ ,  the  minimization  process 
is  greatly  accelerated. 


Minimizing  the  penalty  function  (that  is,  obtaining  the  best  fit) 
can  be  interpreted  geometrically  as  finding  the  bottom  of  a  hypersurface 
whose  dimensionality  is  the  number  of  functional  parameters.  An  initial 
guess  for  the  pk  provides  a  starting  point  for  the  descent.  Taking  an 
initial  guess  that  is  close  to  the  true  bottom  reduces  the  descent 
e  'fort  and  helps  keep  the  descent  path  out  of  undesirable  local  minima 
that  are  far  from  the  true  bottom  (that  is,  the  least  minimum). 
Although  the  computer  calculation  was  performed  in  double  precision,  it 


®R.  Fletcher  and  C.  M.  Reeves ,  Function  Minimization  by  Conjugate 
Gradients,  Computer  Journal,  2.  (^64),  149-154 . 
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was  learned  that  a  minimum  in  the  hypersurface  takes  on  a  rough, 
granular  texture  in  its  shallowest  portions  due  to  roundoff  errors  in 
the  calculation.  This  undesirable  texture  impedes  the  descent  to  the 
desired  minimvsn.  In  general,  the  speed  of  descent  is  slowed  by 

complexity  or  exceptional  shallowness  of  the  hypersurf ace .  Greater 

complexity  is  typically  caused  by  more  parameters  having  higher  poly¬ 
nomial  degree  than  quadratic  {in  a  Taylor  expansion  of  the  penalty 

function)  so  that  the  descent  requires  more  changes  in  direction. 

Shallowness  increases  in  rough  proportion  to  the  number  of  scattered 
data  being  fitted  with  the  functional.  For  example,  if  the  number  of 
data  is  increased  tenfold,  but  the  data  scatter  (due  to  finite  variance) 
is  not  reduced,  the  hypersurface  becomes  roughly  10  times  shallower.  It 
then  becomes  harder  to  find  the  precise  bottom  of  the  shallower 
surface.  In  this  work,  the  penalty  function  hypersurface  is  both 

complex  and  shallow,  and  the  success  of  the  fitting  process  depended 
critically  on  the  efficiency  and  the  robustness  of  the  minimization 

procedures . 

Hie  initial  guesses  for  the  various  sets  of  parameters  p^  were 
usually  taken  from  the  optimal  values  of  pk  for  adjacent  source  neutron 
energy  bands  or  burst  heights.  After  the  optimum  values  for  p^  were 
obtained  for  a  particular  descent,  the  importance  of  each  parameter  was 
tested.  If  elimination  of  a  parameter  caused  a  sizeable  increase  in  the 
penalty  function,  that  parameter  was  deemed  important.  If  a  negligible 
increase  was  observed,  that  parameter  was  deemed  unimportant  and  was 
removed  from  the  functional.  The  functionals  for  several  drivers  were 
considerably  simplified  in  this  way. 

An  example  of  the  effect  of  minimizing  the  penalty  function  is  shown 
for  one  volume  detector  in  table  2.  The  energy  deposition  (or  ioniza¬ 
tion)  due  to  thermal  neutron  capture  in  the  ground  is  given  for  volume 
detector  33,  source  neutron  energy  band  9  (8.19  to  15.0  MeV) ,  and  200-m 
burst  height.  Detector  33  is  600  to  900  m  above  the  ground  at  1200-  to 
1500-m  ground  range  from  the  burst.  Table  2  gives  these  parameters: 

the  time  bin  boundaries, 

fj^,  the  Monte  Carlo  score  in  the  detector, 

f,  the  optimal  global  fit  score  in  the  detector, 

fT,  the  trial  global  fit  score  in  the  detector, 

d,. ,  the  percentage  of  fractional  deviation  ( due  to  the  vari¬ 
ance)  of  the  Monte  Carlo  score  for  this  driver, 

dja ,  the  percentage  of  fractional  deviation  (due  to  the  vari¬ 
ance)  of  the  Monte  Carlo  score  for  the  sun  of  all  drivers. 


r,  the  percentage  fraction  of  the  total  Monte  Carlo 
score,  fij,  due  to  this  driver,  f^  (r  =  100  fij/fij)» 

e,  the  percentage  of  difference  between  the  Monte  Carlo 
score,  f^ j ,  and  the  optimal  fit  score,  f, 

A 

e,  the  percentage  of  "error"  in  the  total  score  contributed 
by  optimal  fit  score  f, 

eT,  the  percentage  of  difference  between  Monte  Carlo  score 
f^j  and  trial  fit  score  fT, 

A  A 

e«p,  the  same  as  e,  but  evaluated  for  the  trial  fit  score 


The  goal  of  the  fitting  process  is  to  make  e  as  small  as  possible  in 
each  time  bin.  When  r  is  larger,  a  better  fit  is  needed  to  make  e  small 
and  thereby  make  e  small*  Generally,  e  is  desired  to  be  roughly  equal 
to  or  smaller  than  d^  for  each  time  bin.  Since  fractional  deviation 
d^ .  (due  to  the  variance)  is  an  estimate  of  the  Monte  Carlo  error, 
difference  e  between  the  fit  and  the  Monte  Carlo  score  should  be  about 
equal  to  d  — .  Thus,  when  d^ ■  is  large,  e  may  be  large;  when  d—  is 
small,  e  should  be  small.  In  ~table  2,  time  bins  11  through  13  haveJ  d^ 
averaging  20  to  25  percent  and  e  averaging  about  10  percent;  time  bins 
15  through  18  have  d^j  averaging  40^  to  50  percent  and  e  averaging  about 
40  percent.  The  values  for  e  and  e  appear  to  fluctuate  approximately 
randomly  about  zero  in  time  bins  8  through  19,  as  would  be  expected  of  a 
good  fit  to  the  data.  The  largest  values  for  eT  and  e  have  been 
reduced  by  the  fitting  process  to  smaller  values  for  e  and  e  • 


Sr*.  " 


TABLE  2.  EFFECT  OF  MINIMIZING  PENALTY 
DETECTOR  33,  NEUTRON  ENERGY 
BURST  HEIGHT 

FUNCTION 
BAND  9  (8. 

FOR  IONIZATION 
,19  TO  15.0  M*V) 

RATE,  VOLUME 
,  200-m 

Tinw 

bln 

Start  of 
tiiM  bln 
(a) 

f 

't 

Oil 

5ti 

r 

• 

• 

•t 

*T 

1 

0.0 

0.0 

0.89020-12 

0.92830-12 

99.9 

9.3 

0.0 

100.0 

0.0 

100.0 

0.0 

2 

0. 10000-06 

0.  70370-08 

0.37490-11 

0.3910D-11 

99.9 

10.8 

0,1 

-99.9 

-0.1 

-99.9 

-0.1 

3 

0. 2150D-06 

0.0 

0.06550-11 

0.90260-11 

99.9 

11.4 

0.0 

100.0 

0.0 

100.0 

0.0 

4 

0.4640D-06 

0.0 

0 1905D-10 

0.19870-10 

99.9 

10.9 

0.0 

100.0 

0.0 

100.0 

0.0 

5 

0.1000D-05 

0.0 

0.40530-10 

0.42280-10 

99.9 

9.9 

0.0 

100.0 

0.0 

100.0 

0.0 

6 

0.21 S0D-05 

0.0 
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-0.4 
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12 
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0.46400-03 
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29.3 

23.7 

80.6 
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-8.7 

-33.5 

-27.0 

14 

0.10000-02 

0.35460-10 

0.  30130-10 

0.1921D-10 

69.2 

58.4 

83.7 

-15.0 

-12.6 

-45.8 

-38.3 

15 

0.21500-02 

0.60820-11 

0.07490-11 

0.6371D-11 

43.6 

36.6 

63.4 

30.5 

19.3 

4.5 

2.9 

16 

0.46400-02 

0.26890-11 

0.381 3D- It 

0.28670-1) 
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28.7 

13.9 

29.5 

4.1 

6.2 

0.9 

17 

0.10000-01 

0.3794D-11 

0 . 16880-11 

0.1283D-11 

39.2 

25.5 

27.7 

-55.5 

-15.4 

-66.2 

-18.3 

18 

0.21500-01 

0.13490-11 

0.7428D-12 

0.5717D-12 

52.4 

15.S 

15.8 

-44.9 

-7.1 

-57.6 

-9.1 

19 

0.46400-01 

0.3319D-12 

0.32620-12 

0.2542D-12 

35.0 

17.5 

7.6 

-1.7 

-0.1 

-23.4 

-1.8 

Motes  t 

t *  t he  ttonta  carlo  score  in  the  detector . 

f  *  eh*  optimal  global  fit  acore  in  the  detector . 

*  the  trial  global  fit  score  in  the  detector • 

dXj  -  th*  percentage  of  fractional  deviation  (due  to  the  variance) 
of  th*  Monte  Carlo  score  for  these  drivers. 

<*X4  *  the  percentage  of  fractional  deviation  (due  to  the  variance ) 
of  Ch*  Monte  Carlo  score  for  the  sum  of  all  drivers. 

r  -  100  fij/fij  . 

•  •  the  percentage  of  difference  between  ty  and  f. 

a  •  re/100. 

eT  *  th*  percentage  of  difference  between  fij  an t  fT. 
eT  •  re j/ 100. 
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6.  BURST  HEIGHT  DEPENDENCE  OF  FIT  PARAMETERS 

Having  obtained  parameter  values  at  four  burst  heights ,  we  desired 
to  develop  smooth  functions  of  burst  height  to  approximate  the  parameter 
values*  Some  scatter  was  evident  in  the  burst  height  dependence.  The 
functions  chosen  should  avoid  reproducing  the  scatter*  but  capture  the 
important  trend  of  the  data.  TVo  burst  height  functions  were 

selected:  (1)  the  constant  function,  independent  of  burst  height,  and 

(2)  the  rational  function 

Ci  +  C2h  +  C3h2 

g(h)  **  -  • 

1  +  C4h2 

where  h  is  the  burst  height  and  Ci  are  fit  coefficients.  By  inspection 
of  the  parameter  values  at  the  four  burst  heights,  the  author  determined 
which  of  the  two  burst  height  functions  was  appropriate  for  a 
parameter.  Many  parameters  were  clearly  independent  of  burst  height. 
Other  parameters  displayed  such  severe  scatter  that  it  was  necessary  to 
use  merely  an  average  value  for  all  burst  heights.  On  the  other  hand, 
many  parameters  showed  well-defined  trends  and  were  fitted  with  rational 
functions . 

Fitting  the  four  values  of  a  parameter  with  a  rational  function  was 
not  a  trivial  problem.  In  its  application  to  real  EMP  prediction 
problems,  the  rational  function  would  be  evaluated  at  other  burst 
heights.  Therefore,  the  behavior  of  the  rational  function  must  be 
carefully  controlled  for  all  possible  burst  heights.  FOr  example, 
simply  fitting  the  four  parameter  values  with  rational  function  g(h) 
would  reproduce  the  data  exactly,  but  might  introduce  poles,  crossovers, 
and  similar  erratic  behavior  at  other  burst  heights.  To  avoid  this 
possibility,  several  strategems  were  adopted.  First,  the  fit  was  to  be 
done  by  use  of  Courant's  well-known  penalty  function  method  for  optimiz¬ 
ation  with  nonlinear  constraints,  and  the  nonlinear  constraints  would  be 
defined  so  as  to  inhibit  the  erratic  behavior.  Second,  the  rational 
function  would  be  made  to  fit  not  the  four  data,  but  the  piecewise 
linear  function  p  interpolating  the  four  data. 

The  penalty  function  used  was 

Q(ci?r)  -  £  [g^)  -  p^)]  +  r(N!  ♦  N2  ♦  N3)  . 


The  hj  are  about  40  burst  heights  on  the  interval  from  1  to  300  m. 

,  n2,  and  Nj  are  positive  semidefinite  functions  that  are  zero  when  a 
corresponding  constraint  i3  satisfied.  In  the  optimization  process,  Q 
is  minimized  for  successively  larger  values  of  r  by  using  as  each  start¬ 
ing  point  the  previous  optimum  values  for  .  As  r  is  increased  to 
infinity,  Q  becomes  optimal  for  the  constraints  represented  by  Ni,  N2, 
and  Nj .  The  specific  nonlinear  constraints  were  these: 

a.  For  complex  h,  no  poles  of  g(h)  in  the  right  half-plane  of  h, 

NJ  =  108[min(o,C4)]2  , 

b.  No  crossover  to  negative  values  of  g(h)  for  large  h, 

N2  =  108[min(0,C3)]2/p(200  m)  , 

c.  No  negative  values  for  g(h)  for  small  h, 

N3  =  [min(o,Ci) ]2/p{200  m)  . 

The  minimization  of  Q  was  done  for  each  value  of  r  by  adjusting  the  four 
C^.  The  descent  algorithms  described  in  section  5  were  used.  For  burst 
heights  above  200  m,  most  driver  parameters  tend  toward  limiting  values 
not  much  different  from  the  value  at  the  200-m  burst  height.  An  excep¬ 
tion  is  the  amplitude  parameter  for  a  ground- induced  source.  This 
parameter  should  decay  monotonically  to  zero  as  the  burst  height  in¬ 
creases.  For  burst  heights  above  250  m,  the  three  high-energy  ground 
driver  amplitudes  are  made  to  decay  by  a  factor  of  two  for  each  addi¬ 
tional  100  m.  Similarly,  the  three  low-energy  ground  driver  amplitudes 
are  made  to  decay  by  a  factor  of  1.2  to  4,  depending  on  the  source 
neutron  energy,  fcr  each  additional  100-m  increase  in  burst  height  over 
250  m. 

7.  COMPARISON  WITH  LIQUID  AIR  EXPERIMENT 

Sidhu  et  al9  measured  neutron  and  secondary  gamma  ray  transport 
through  a  liquid  air  sphere  for  a  14-MeV  neutron  source.  Their  129.3-cm 
sphere  radius  was  equivalent  to  an  atmospheric  range  of  850  m.  They 
published  (their  fig.  7)  for  secondary  gamma  rays  a  value  of  about  1.3  x 
10"10  cm2 -rad  (tissue) /source  neutron  for  the  4wr2  dose  at  850-m  equiva¬ 
lent  range  in  air.  The  corresponding  value  for  the  NEMP  code  energy 
band  9  (8.19  to  15.0  MeV)  high-energy  air  driver  is  1.5  x  10~10  cm2-rad 
(air).  One  rad  in  tissue  and  one  rad  in  air  differ  by  only  about  1 
percent,  so  these  gamma  doses  differ  by  about  15  percent. 

9G.  5.  Sidh u,  W,  E.  Farley,  L.  F.  Hansen,  T.  Komoto,  B.  Pohl,  and  C. 
Wong,  Transport  of  neutron  and  Secondary  Gamma  Radiations  Through  a 
Liquid  Air  Sphere  Surrounding  a  14-MeV  Neutron  Source,  Nuclear  Science 
and  Engineering,  66_  ( June  1978),  428-433, 
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The  NEMP  neutron  dose  in  air  cannot  be  compared  easily  with  the 
measured  neutron  tissue  dose  because  neutron  kerma*  factors  are  very 
energy  dependent  in  the  range  from  1  to  14  MeV.  The  neutron  energy 
spectrum  at  the  850-m  observer  for  the  NEMP  driver  is  unavailable  to 
allow  a  conversion  fron  neutron  air  dose  to  neutron  tissue  dose. 

Shown  in  figure  7  of  Sidhu  et  al9  are  results  of  TARTNP  computer 
code  calculations  of  the  gamma  tissue  dose.  The  TARTNP  code  is  a  dis¬ 
crete  ordinates  radiation  transport  code  used  at  the  Lawrence  Livermore 
Laboratories.  The  results  show  the  4irr2  dose  to  a  2800-m  range. 
Differences  between  the  NEMP  driver  and  TARTNP  results  appear  to  be 
small  (a  few  percent)  from  850  to  2500  m.  The  NEMP  code  results  are 
roughly  15  pecent  lower  at  the  500-m  range . 

8.  COMPARISON  WITH  EARLIER  APPROXIMATION  FOR  IONIZATION  BY  14-MeV 

NEUTRONS 

Longley  et  al5 * * * 9  developed  approximations  for  BMP  drivers  for  use  in 
their  LIMP  2  surface  burst  BMP  environment  prediction  code.  They  used 
Longmire's  analytic  first-scatter  theory  and  many-scatter  diffusion 
theory  to  describe  14-MeV  neutron  direct  ionization  through  elastic 
scattering  and  (n,p)  reactions.  Comparison  of  their  14-MeV  neutron 
driver  with  the  corresponding  NEMP  driver  (energy  band  9,  direct  neutron 
ionization)  showed  rough  agreement  and  some  differences.  For  dose  in 
air,  the  NEMP  results  exhibited  a  10-percent  longer  attenuation  length 
at  all  ranges  and  a  roughly  50-percent  greater  amplitude  (factoring  out 
the  range  dependence).  Because  detailed  reaction  cross  sections  were 
not  used  for  this  LEMP  2  driver  (per  C.  L.  Longmire,  Mission  Research 
Corp.),  some  disagreement  is  expected. 

Since  the  LEMP  2  drivers  are  independent  of  the  source  neutron 
energy  spectrum  (except  for  14-MeV  neutrons),  a  comparison  with  the  nine 
energy  band  NEMP  drivers  is  a  complicated  undertaking  beyond  the  scope 
of  this  report. 


SH.  J.  longley,  C.  L.  longmire,  J.  S.  Malik,  R.  M.  Hamilton,  R.  N. 

Marks,  and  K.  S.  Smith,  Development  and  Testing  of  LEMP  2,  A  Surface 

Burst  EMP  Oode  (V),  Mission  Research  Oorp.,  Santa  Barbara,  CA,  DNA  4097T 

(December  1976),  (CONFIDENTIAL) 

9G.  S,  Sidhu,  W.  E,  Farley,  L,  F.  Hansen,  T,  Komoto,  B,  Pohl ,  and  C. 
Wong,  Transport  of  Neutron  and  Secondary  Gamma  Radiations  Through  a 
Liquid  Air  Sphere  Surrounding  a  14-Mev  Neutron  source.  Nuclear  Science 
and  Engineering,  66_  (June  1978),  428-433. 

*The  total  kinetic  energy  of  directly  ionizing  particles  ejected  by  the 
action  of  indirectly  ionizing  radiation  per  unit  mass  of  specified 
material . 


9.  SELECTED  RESULTS 

Time  histories  of  some  NEMP  neutron- induced  EMP  drivers  are  por¬ 
trayed  in  figures  2  to  13.  Radial  Compton  current,  theta  Compton 
current,  and  ionization  rate  are  plotted.  Results  are  shown  for  two 
observers  0  and  500  m  above  the  ground  for  each  of  two  burst  heights  of 
1  and  200  m.  The  ground  range  from  the  burst  to  the  observer  is  1000 
m.  All  results  are  for  a  source  of  one  neutron  in  the  energy  band  of 
8.19  to  15.0  MeV.  Although  drivers  for  other  neutron  energy  bands  have 
different  amplitudes,  attenuation  lengths,  burst  height  dependence,  air- 
ground  interface  effects,  etc.,  these  results  are  typical. 

In  the  figures,  each  driver  is  represented  by  a  broken  line  with 
identifying  symbols.  The  "F"  denotes  the  high-energy  ground  driver 
(ground  inelastics);  the  "I"  denotes  the  high-energy  air  driver  (air 
inelastics);  the  "G"  denotes  the  low-energy  ground  driver  (ground  cap¬ 
tures);  the  "A"  denotes  the  low-energy  air  driver  (air  captures);  and 
the  "N”  denotes  neutron  direct  ionization  including  charged  particle 
ionization.  The  total  for  all  drivers  is  represented  by  a  solid  line 
with  identifying  "T"  symbols.  Each  curve  is  the  magnitude,  or  absolute 
value,  of  the  driver  or  total. 

Several  general  observations  can  be  made  concerning  burst  height  and 
observer  height  variations.  High-energy  ground  drivers  are  signifi¬ 
cantly  weaker  for  greater  burst  heights.  Low-energy  ground  drivers 
decay  more  slowly  with  increasing  burst  height  than  do  high-energy 
ground  drivers  because  source  neutrons  penetrate  to  greater  ranges  while 
decaying  to  lower  energies.  Low-  and  high-energy  air  drivers  and  the 
neutron  direct  ionization  driver  all  are  roughly  twice  as  strong  for  the 
200-m  burst  height  as  for  the  In  burst  height  because  roughly  half  of 
the  source  neutrons  are  absorbed  by  the  ground  for  the  lower  burst 
height. 

The  drivers  for  theta  Compton  current  are  all  stronger  for  observers 
closer  to  the  ground  surface  because  of  the  asymmetry  presented  by  the 
air-ground  interface.  High-energy  and  low-energy  air  drivers  for  radial 
Compton  current  and  ionization  rate  are  relatively  insensitive  to 
observer  height.  The  time  dependence  of  the  neutron  direct  ionization 
driver  is  different  for  the  two  observer  heights  considered.  This 
driver  is  somewhat  shorter  in  duration  near  the  ground  because  many 
later-arriving  multiply  scattered  neutrons  tend  to  become  trapped  in  the 
ground  and  do  not  contribute  to  the  driver  in  the  air.  The  low-energy 


ground  driver  for  ionization  rate  is  somewhat  stronger  for  the  higher 
observer,  more  noticeably  so  for  the  200-m  burst  height  than  for  the  1-m 
burst  height.  They  differ  because  the  secondary  gamma  rays  (produced  by 
neutron  capture  in  the  ground)  that  drive  the  ionization  rate  at  the 
observer  must  penetrate  more  soil  to  reach  the  shallower  observer.  The 
low-energy  ground  driver  for  radial  Compton  current  is  even  more 
strongly  decreased  for  observers  near  the  ground  because  the  net  Compton 
current  flow  follows  the  direction  of  the  gamma  ray  flux,  which  is 
predominantly  up  out  of  or  down  into  the  ground  (perpendicular  to  the 
radial  direction).  Higher  off  the  ground,  more  gamma  rays  with  a  nearly 
radial  direction  arrive  from  ground  regions  close  to  the  burst. 

All  drivers  contributing  to  the  radial  Oompton  current  have  the  same 
sign  (negative).  The  theta  Compton  current  drivers  may  have  different 
signs  or  even  change  sign.  (By  convention,  positive  is  downward.)  The 
theta  Compton  currents  for  high-energy  ground  and  low-energy  air  drivers 
are  positive  and  negative,  respectively.  The  theta  Oompton  current  for 
the  low-energy  ground  driver  is  generally  negative  at  this  range  from 
the  burst  because  of  a  secondary  gamma  ray  "fountain"  effect.  This 
effect  is  due  to  many  secondary  gamma  rays  rising  upward  from  the  ground 
close  to  the  burst  and  being  turned  down  toward  the  ground  after  a  few 
Compton  collisions  in  the  air.  However,  for  observers  close  to  the 
ground  and  close  to  the  burst,  enough  upward- going  secondary  gamma  rays 
are  produced  beneath  the  observer  to  reverse  the  theta  Compton  current 
from  negative  to  positive  for  a  period  of  time  (fig.  9).  The  theta 
Compton  current  for  high-energy  air  reactions  was  set  identically  to 
zero  because  it  was  about  two  orders  of  magnitude  smaller  than  the 
corresponding  radial  Compton  current  and  because  the  available  Monte 
Carlo  statistics  for  those  data  were  so  poor  that  a  reliable  fit  was 
unlikely. 

The  time  history  of  the  high-energy  ground  drivers  for  theta  Compton 
current  and  ionization  rate  deserves  some  explanation.  It  may  consist 
of  two  pulses:  one  from  the  initial  arrival  of  neutrons  at  the  ground 
beneath  the  burst  and  the  other  from  the  arrival  of  neutrons  at  an  area 
roughly  midway  between  the  burst  and  the  observer.  The  secondary  gamma 
rays  produced  in  the  ground  have  an  attenuation  length  considerably 
greater  than  the  mean  free  path  of  the  neutrons.  This  greater  length 
causes  contributions  from  reactions  (per  unit  ground  surface  area) 
beneath  the  burst  to  be  more  important  than  contributions  from  reactions 
very  near  the  observer.  However,  gamma  rays  originating  in  the  ground 
closer  to  the  observer  are  less  attenuated  by  the  ground  because  of 
their  higher  angle  of  departure  from  the  ground.  This  complicated 
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behavior  was  approximated  by  the  two  pulses  mentioned.  The  approxi¬ 
mation  gives  reasonable  agreement  with  the  Monte  Carlo  data,  but  too 
much  importance  should  not  be  attached  to  the  detailed  structure  of  the 
time  history  approximation,  expecially  since  the  data  time  bins  were 
one-third  of  an  order  of  magnitude  wide.  The  second  pulse  of  this  two- 
pulse  structure  is  not  observed  in  the  data  (and  the  driver)  for  radial 
Compton  current.  Also,  the  second  pulse  is  not  observed  in  the  data 
(and  the  driver)  for  ionization  rate  if  the  burst  height  is  very  near 
the  surface  (fig.  4  and  7)  or  if  the  observer  is  very  high  (fig.  13). 


Figure  2.  Radial  Compton  current  for  1-m  burst  height,  0-m 
observer  height,  and  1000-m  ground  range. 
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Figure  3.  Theta  Compton  current  for  1-m  burst  height,  0-m 
observer  height,  and  1000-m  ground  range. 
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Figure  8.  Radial  Compton  current  for  200-a  buret  height,  0-m 
observer  height,  and  1000-m  ground  range. 
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height,  and  1000-n  ground  range. 
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Figure  13.  Ionization  rate  for  200-m  burst  height,  500-m  observer 
height,  and  1000-m  ground  range. 


10.  CONCLUSIONS 


The  NEMP  neutron- induced  EMP  driver  package  represents  the  state  of 
the  art  for  predicting  EMP  drivers  from  near-surface  bursts,  allowing 
flexible  representation  of  any  weapon's  neutron  energy  spectrum  and 
choice  of  any  burst  height  in  the  near-surface  regime.  The  package 
consists  of  separate  space-  and  time- dependent  drivers  for  radial  and 
theta  Compton  current  and  ionization  rate,  for  high-  and  low-energy 
reactions  in  the  air  and  the  ground,  for  nine  source  neutron  energy 
bands,  for  any  near-surface  burst  height. 

This  report  has  described  the  development  of  this  EMP  driver  package 
for  neutron- induced  sources,  for  use  in  the  NEMP  code  for  EMP  environ¬ 
ment  prediction  for  near-surface  nuclear  bursts.  The  development  was 
based  on  the  concept  of  obtaining  a  globed  fit  of  certain  well-chosen 
functionals  to  Monte  Carlo  transport  data  over  the  entire  set  of  volume 
detectors  and  time  bins.  The  global  fit  obtained  in  this  manner  offered 
several  advantages  over  fitting  with  polynomials,  including  the 
effective  elimination  of  Monte  Carlo  "noise"  and  the  possibility  of 
reduction  of  error  below  the  level  of  Monte  Carlo  statistical  error. 
The  fit  also  involved  far  fewer  degrees  of  freedom.  The  global  fit  was 
obtained  by  minimization  of  a  penalty  function  representing  the  badness 
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of  fit.  Because  of  the  large  amount  of  computer  time  involved,  highly 
efficient  minimization  algorithms  were  employed,  the  difficulty  of  the 
minimization  depended,  among  other  things,  on  the  variance  of  the  Monte 
Carlo  data,  which  averaged  about  20  to  30  percent.  Ifce  burst  height 
dependence  of  the  fit  parameters  was  explicitly  obtained  for  burst 
heights  from  1  to  250  m.  Extrapolation  methods  were  devised  for  burst 
heights  greater  than  250  m. 

Comparison  of  ionization  from  gamma  rays  induced  by  a  14-MeV  neutron 
source  with  experimental  measurements  gave  good  agreement.  Comparison 
with  the  LEMP  2  direct  neutron  ionization  driver  for  14-MeV  neutrons 
revealed  significant  differences.  Typical  results  were  discussed  for 
two  burst  heights  and  two  observer  heights,  for  the  8. 19-  to  15.0-MeV 
source  neutron  energy  band  drivers. 
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APPENDIX  A.— A  TWO-POINT  QUADRATURE  SCHEME 


Integrating  a  driver  functional  over  the  entire  volume  of  a  volume 
detector  can  be  very  time-consuming  if  unsophisticated  quadrature 
methods  are  used.  lhe  volume  detector  concerned  is  a  square  toroid, 
that  is,  a  toroid  with  a  square  cross  section  in  coordinates  (x,z,4). 
The  toroid  is  bounded  by  X^.X-.Z^,  and  Z  ,  where  x  is  the  ground  range 
from  the  coordinate  center  ana  z  is  the  neight  above  the  ground.  The 
driver  has  azimuthal  symmetry  around  the  center  of  the  toroid  (x  *  0)  so 
that  an  integration  in  independent  variable  $  is  trivial.  Thus,  an 
integral  over  detector  volume 


r  r  f« 

Jxi  Jzi  -Jo 


x,z)  (x 


9  O \ l/2 
2  +  z2)  dx  dz  d$ 


becomes 


I  =  2ir 


itT  +  *2)1/2  d*  dz  . 


If  f(x,z)  describes  a  driver  functional,  the  salient  spatial  gradients 
of  f ( x , z )  are  (1)  in  the  vertical  (z  >  0)  direction  due  to  proximity  of 
the  ground  and  (2)  in  the  radial  direction  from  the  burst  due  to 
geometric  and  atmospheric  attenuation.  A  particularly  accurate  Gaussian 
quadrature  rule  would  use  the  geometric  and  atmospheric  attenuation  as  a 
weight  function. 

Our  goal  is  a  Gaussian  quadrature  rule  for  the  integral 

Tx2  fz2 

J  =  I  I  v(x,z)w(x,z;a)  dx  dz  , 

JZl 

where  v  is  a  function  of  x  and  y  and  w  is  a  weight  function  defined  by 


w(x,z;a)  =  exp(-ar)/(l  +  r2)  , 

r  *  [x2  +(z  -  h)2]1^2  . 


The  burst  height  is  h.  The  advantage  of  this  weight  function  rests  on 
the  fact  that  it  represents  the  actual  spatial  dependence  ( for  a 
suitable  choice  of  a)  better  than  a  low-order  power  series  would. 
Parameter  o  is  essentially  the  atmospheric  attenuation  coefficient  (the 
reciprocal  of  the  attenuation  length)  and  is  available  for  use  during 
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the  calculation  of  a  driver  integral*  Thus,  we  may  exploit  our  know¬ 
ledge  of  a  to  construct  the  quadrature  based  on  weight  function 
w(x,z;a).  In  computer  program  BIGFIT,  we  actually  first  construct  a 
family  of  quadrature  rules  for  a  range  of  values  a, *  Then  when  a  driver 
integral  is  being  calculated,  the  nearest  to  a  is  determined,  and  its 
corresponding  quadrature  rule  is  used.  In  this  way,  the  construction  of 
the  quadrature  rules  is  kept  out  of  the  inner  program  loops. 

For  the  J  integral  quadrature  rule,  we  seek  quadrature  points 
(XjfZj)  and  (Xj,z  )  with  weight  coefficients  w^  and  w^,  respectively. 
We  require  exact  Integration  of  Taylor  series  terms  ror  x®  (and  z®), 
xl,z*,z2,  and  z2.  This  requirement  yields  the  following  conditions  on 

WV*!'  and  W2S 


wi  +  W2  =  Mqo  , 
wlxl  +  w2*l  *  M10  * 


wlzl  +  W2Z2  =  M01  » 
wlz2  +  W2z2  »  MO 2  # 
wiz2  +  W2z2  -  MO 3  , 


where 


z;a)  dx  dz 


This  quadrature  is  separable  in  x  and  z  and  easy  to  obtain.  Obviously, 


X1  =  M10/M00* 


Define  the  polynomial 


P  **  (z  -Zj)(z  “Z^)  ■  0  , 


or 


P*z2-Az+B-0. 
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Multiply  the  preceding  equations  by  A  and  B  and  add 
wjB  +  W2B  —  BM00 
W}(Azi)  +  W2(AZ2)  “  AM01 
Wi(z2)  +  »2(z|)  “  m02 

Wi<P)  +  W2(P)  “  AM01  +  BMqo  +  M02 

w^z^(B)  +  W2Z2(B)  *  BMqj 
wjZj (AZ} )  +  W2Z2(AZ2>  “  AM02 
WiZ^z^)  +  '<l2zz(^')m  m03 

w^tP)  +  w2z2(P)  *  AMq2  +  BMqi  +  »03  * 

Since  P  *  0,  we  have 

AM0i  +  BMqo  +  M02  “  0  * 

AM 02  +  BMqj  +  Mq3  ■  0  , 

or  solving  for  A  and  B, 

A  *  (M03M00  ”  Mq2m0i)/c  ' 

B  -  (Mj2  -  M01M03)/C  , 

where 

C  -  -  MO0M02  • 


Therefore, 


.j  -  [-A  -  (A2  -  4b)1/2]/2  , 

-  [-*  ♦  (A*  -  4B)l/2]/2 
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.Ml  ,'BP 


Then 

W1  +  w2  "  M00  • 
wlzl  +  W2Z2  *  *01  ' 

or 

wi  *  (Mqi  -  Moo*2)/(z1  “  z2)  * 

w2  *  Mqq  -  Wj  . 


Suppose  we  seek  to  Integrate 

fz2 

,z)  dx  dz 

where  y(x,z)  has  dominant  spatial  variation  w(x,z;a).  He  rewrite  K  as 
K  *  //  [)f^x«ziajl h(x,z;o)  dx  dz 


nz2 

U<x,j 

. 


and  apply  the  quadrature  rule  developed  to  get 

H*  (X1  *zl )  1  r *2)  ] 

*  “  Wl l-H^x  1 »z  j  ju]  -1  +  W2Lw(xi,z2,a)  J 

It  is  best  to  combine  w(x  ,z  ;a)  with  wi  and  w(x,,z  »<*)  with  w  to 
obtain  the  modified  rule 

*  ’  V(Vh)  *  V(V*2)  • 


where 

A 

W1  “  W1^W(X1#*1,°)  ' 

A 

w2  -  Wg/w^Xj  »z2 ;a^  . 

A  A 

Thus,  a  unique  quadrature  rule  (x^z^.z  ,Wj,  and  w^)  was  obtained 
for  each  combination  of  volume  detector  and  before  commencing  the 
descent  algorithm  to  fit  a  driver  functional  to  Monte  Carlo  data.  Then 
during  the  descent,  the  numerical  integration  of  K-like  integrals 
proceeded  extremely  rapidly. 
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APPENDIX  B. — DESCENT  ALGORITHMS  USED 

This  appendix  lists  the  FORTRAN  version  of  the  descent  algorithms 
used  to  search  for  a  minimum  of  a  function* 
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SUBROUTINE  TALC6(FN»X. TOL. B.C.XOLD. STEP. VAL.M.NFN) 

C 

C  —  BY  KILLIAN  T.  HYATT*  JR.,  19TB. 

C 

C  USING  THE  NETHOD  OF  CONJUGATE  GRADIENTS*  THIS  SUBROUTINE  SEARCHES  FOR 
C  THE  NININUN  OF  THE  SCALAR  FUNCTION  FN  *  GIVEN  A  STARTING  POINT  X  . 

C  X  IS  A  VECTOR  OF  LENGTH  N  .  AN  ELENENT  OF  X  IS  NOT  ALTERED 
C  DURING  THE  SEARCH  IF  THE  CORRESPONDING  ELENENT  OF  TOL  IS  NON- 
C  POSITIVE.  OTHERHISE*  THE  SEARCH  CONTINUES  UNTIL  THE  FUNCTION  DOES 
C  NOT  OECREASE  HHEN  X  IS  CHANCED  BY  AS  NUCH  AS  TOL«STEP  *  NHERE  THE 
C  LARGEST  ELENENT  IN  STEP  IS  UNITY. 

C 

C  THE  SEARCH  NAY  BE  ABORTED  AFTER  ABOUT  NFN  EVALUATIONS  OF  THE 
C  FUNCTION  FN  IF  NFN>0  .  OR  AFTER  -NFN  GRADIENT  EVALUATIONS  IF 
C  NFN<0  • 

C 

C  THE  FUNCTION  FN  HUST  BE  DECLARED  EXTERNAL  IN  THE  CALLING  PROGRAM. 
C  X  .  TOL  .  N  .  AND  NFN  MUST  BE  DEFINED  HHEN  THIS  SUBROUTINE  IS 
C  CALLED.  THE  SUBROUTINE  RETURNS  THE  MINIMUM  LOCATION  IN  X  *  THE 
C  MINIMUM  FUNCTION  VALUE  IN  VAL  .  AND  THE  NUMBER  OF  FUNCTION 
C  EVALUATIONS  IN  ABSINFN I  .  NFN  IS  RETURNED  NEGATIVE  IF  THE  SEARCH 
C  IS  ABORTED  AS  DESCRIBED  ABOVE.  THE  ARRAYS  B  *  G  *  XOLO  *  ANO 
C  STEP  OF  LENGTH  M  ARE  USED  FOR  INTERMEDIATE  RESULTS.  TOL  AND  M 
C  ARE  RETURNED  UNCHANGED. 

C 

C  THE  PROGRAM  CALLS  THE  SUBROUTINE  QUAD  TO  LOCATE  A  MINIMUM  ALONG  A 
C  GIVEN  DIRECTION. 

C 

C  AN  ELEMENT  OF  X  MAY  BE  COMPELLED  TO  REMAIN  POSITIVE  DURING  THE 
C  SEARCH.  BY  USING  THE  TNO  STATEMENTS  (MARKED  KITH  -CCCCC"!  IN  THE 
C  PROGRAM  BELOMX 

C  TOLtI l*-TOLf I) 

C  GO  TO  5 

c 

C  THE  METHOD  ESTIMATES  THE  FUNCTION  GRADIENT  BY  FINITE  DIFFERENCES  OVER 
C  AN  INTERVAL  TQL*SLASH  .  SLASH  HAS  BEEN  SET  TO  O.OI  IN  A  DATA 
C  STATEMENT.  ROUNOOFF  ERRORS  HILL  CAUSE  POOR  ESTIMATES  FOR  THE 
C  GRADIENT  IF  TOL  IS  TOO  SMALL  OR  IF  MACHINE  PRECISION  IS 
C  INSUFFICIENT.  (DOUBLE  PRECISION  IS  RECOMMENOED  FOR  COMPUTERS  KITH 
C  DEFAULT  SHORT  PRECISION  MORDS.I 
C 
C 

DIMENSION  XIM I .TOL (M I «BIMI *G(M I .XOLDIM I .STEP IM) 

DATA  SLASH  / 0.01/ 

C 

C  INITIALIZE  GRADIENT  AND  DIRECTION  ARRAYS  AND  OTHER  QUANTITIES. 
NAXFN«NFN 
NFN«0 
IPASS-0 

S  CONTINUE 

DO  10  1*1 »N 
1(11*0.0 
10  G( I) *1 .0 

VAL*FN(XI 
NFN*NFN«1 
C 
C 

C  MINIMIZATION  LOOP  BEGINS  HERE.  CHECK  FOR  MAX  FN  EVALUATIONS. 

100  CONTINUE 
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IPASS*1PASS«1 

1FI1PASS.OT.-MAXFN.ANO.MAXFN.LT.O)  CO  TO  900 
1FINFN.CT.MAXFN.AND.MAXFN.0T.0)  00  TO  900 
C 

C  FUST.  MOVE  X  TO  XOLD. 

SUM0LD*0.0 
00  110  1*1, M 
XOLD I I 1*XI II 

110  SUHOLO*SUMOLO«Ctll«*2 

1FISUN0LD.EQ.0.0)  00  TO  S 
FNOLD*VAL 

C  NEXT,  SET  OAAOIENT  ANO  MEM  SEARCH  DIRECTION. 

SUM*0.0 
DO  120  1*1 .N 
0111*0.0 

lFITOLIU.LE.O.Ol  OQ  TO  120 
Kill *XOLD ( I l+SL ASH*TOL I I 1 
VAL*FNIXI 
NFN*NFN*1 
XIII *XOLO( 11 
01 1 1  * IVAL-FNOLO 1/SLASH 
120  SUM*SUM*Ct 11**2 
SUM-SUH/SUMOLD 
DO  130  1*1. N 
Rill *“011 1 +SUM*BI 1 1 
130  STEP  II 1*BI 11 
C 

C  SCALE  STEP  SO  EACH  ELEMENT  1$  .LE.  1.0. 

SAVE *0.0 
DO  140  1*1. N 
TEMP*ABSISTEPI1 1 1 
IF  IT EHP.CT .SAVE )  SAVE*TEMP 
140  CONTINUE 

1FISAVE.NE.0.01  SAVE *1 .0 /SAVE 
C 

C  TARE  TRIAL  STEP. 

DO  ISO  1*1. N 

STEP  II 1*STEPI 1 l*TOLI I 1 •SAVE 
C  LOCKUP  PARAMETER  TENDINO  TO  ZERO. 

IF (STEP (I  I .OE .0 .0 .OR .XDL Dill .OT .T0LII11  00  TO  ISO 
CCCCC  TOLU  1*— TOLI 1 1 
CCCCC  00  TO  S 
ISO  X|1I*XOLOI11«STEPI11 
VAL*FNIXI 
NFN*NFN«1 

IFIVAL.OE .FN0LD1  00  TO  910 
C  MINIMUM  IMPROVED.  00  SEARCH  FOR  MINIMUM. 

CALL  QUADIFM.M.OLOVAL.VAL.XOLD.STEP.X.NFNI 
CCCCC  PRINT  9T.VAL 
97  FORMAT I  *  MEN  VAL**.f20.B) 

CO  TO  100 

c 

c 

C  FINISHED. 

900  NFN—NFN 

DO  905  I  * l.M 
90S  T0LI1I*ABSIT0L«III 
RETURN 

910  CONTINUE 

IF IIPASS.OT.S)  00  TO  91S 
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C  TRY  THROWING  OUT  A  BAO  PARAMETER  IF  IN  EARLY  STASES  OF  NIRINI2AT10N. 
DO  914  1*1  «M 

IFISTEPI1  )  .EQ.O.Ol  SO  TO  914 

SAVE -STEP! II 

STEPI1I-0.0 

XI1)*X0LDI1) 

VAL-FNIX1 

NFN*NFN»1 

IFIVAL.SE. FMDLD1  SO  TO  913 

CALL  QUAD  I FN ,H , DLDVAL , VAL .XDLD.STEP.X tNFNI 

•111*0.0 

Slll-0.0 

SO  TO  100 

913  STEP (1 I *$AVE 

Kill *XOLD( 1I*STEP( 1} 

914  CONTINUE 

915  CONTINUE 

C  RESTORE  X  AND  VAL. 

VAL*FNOLD 
DO  920  1*1. N 
TOLIll*AiS(TOLCIII 
920  XI II *XOLD 1 11 

RETURN 
ENO 


APPENDIX  B 

SUBROUTINE  QUAD IFN.MPARAM.FNOLD, VAL, XOLD .STEP, X.NFN) 

—  BY  KILLIAN  T.  HYATT.  JR.,  19TB. 

THIS  SUBROUTINE  FINOS  THE  NIM1NUN  OF  THE  FUNCTION  FN  IN  THE  DIRECTION 
STEF,  LIVEN  FNOLD>FNIXOLDI,  VAL»FNIX),  X«XOLD«STEP,  AND  RPARAN  ITHE 
LEN6TH  OF  VECTORS  X,  XOLD,  AND  STEFI.  A  NECESSARY  CONDITION  IS  THAT 
VAL  <  FNOLD.  NFN,  A  COUNTER  FOR  FUNCTION  CALCS.  IS  INCRERERTED  FOR 
EACH  EVALUATION  OF  FN.  THE  SUBROUTINE  RETURNS  A  HINIRUN  VAL  AND 
AR6UNENT  X  SUCH  THAT  VAL-FNtXI,  AND  STEF  SUCH  THAT  STEP-X-XOLO. 

THE  AL60R1THN  USES  QUADRATIC  INTERPOLATION  ACCELERATED  BY  A  LDtARlTM- 
NIC  BISECTION-TYPE  SEARCH. 

D1NENS10N  XOLD INF ARAN) ,STEPINPARAN),X INF ARAN  I 
I  COUNT  UNITS  THE  NUNBER  OF  FUNCTION  EVALUATIONS  TO  A  RUNNER 
CONVENIENT  FOR  THIS  APPLICATION.  THE  USER  NAY  CHANCE  IT. 

1 COUNT *AMAX1 I3.0,SQRT(FL0ATINPARANI) I 

IC»0 

SA«0  . 

SB*1 .0 

FA-FNOLD 

FB-VAL 

LOOK  FOR  A  FC.6E.FB. 

SC«S8*10.0 

DO  20  l-l.HPARAM 

XH)*X0LD(I)«STEPI1I*SC 

CONTINUE 

FC'FNIX) 

NFN>NFN»1 

IFIFC.6E.FB)  60  TO  l&O 
TRY  A6A1N. 

FA»FB 
FB«FC 
SA-SR 
SB-SC 
60  TO  10 

BE6IN  MODIFIED  QUADRATIC  INTERPOLATION  FOR  NIRIRUN. 

D  CORT INUE 
IC«!C«I 

IFI1C.6T.I COUNT )  60  TO  170 
TENP  *(SC— SBI/ISB-SA) 

IFITENP.6T .6.0)  60  TO  110 
1FITENP.LT.0.2S)  60  TO  120 
IFIFA.EQ.FB.OR.FB.EQ.FC)  60  TO  170 
CA«IIFA-f B)/ISA-SBI— IFB-FCI/ISB— SC ll/ISA-SC) 

IF ICA.EQ.O .0 )  60  TO  110 
CD» IFA-FBI /ISA-SBI-CA* (SA«SB ) 

$X>— CD/12. 0*CAt 
60  TO  130 
)  CONT INUE 

EITHER  C0NVER6IN6  VERY  HELL,  OR  POORLY.  USE  L06ARITHNIC  SEARCH. 

TENP -SORT  I TENP 1-1.0 
SX>SB«(SB-SA)*TEMP 
60  TO  130 

)  TENP "SORT  1 1.0/TEMP  I— 1.0 
SX  »  SR-I SC -SB ) *T  ENP 
)  CONTINUE 

IF ISX.EQ.SB.OR.SX.LE .SA.OR.SX .6E.SC)  60  TO  170 
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C  EVALUATE  FUNCTION  AT  INTERPOLATED  POINT. 
DO  135  1*1, HP At AN 
XIII  *XOLD I IMSTEPI 1 1  *SX 
135  CONTINUE 
FX*FNIXI 
NFN*NFN«1 

IFIFX.LT. FBI  CO  TO  150 
C  NOT  NEM  NINIHUN. 

IF ISX.CT .SB)  CO  TO  ICO 
SA-SX 
FA«FX 
CO  TO  100 
ICO  sc*sx 
FC*FX 
CO  TO  100 
ISO  CONTINUE 
C  IS  NEN  NIN1NUM. 

IFISX.LT. SB)  CO  TO  ICO 
SA*SB 
FA*FB 
SB*SX 
FB*FX 
CO  TO  100 
ICO  SC«SB 
FC*FB 
SB*SX 
FB-FX 
CO  TO  100 
170  CONTINUE 
VAL«FB 

00  ISO  1*1 (NPARAM 
STEP! 1 1 “STEP  II MSB 
XII) “XOLD I II ♦STEP 111 
1B0  CONTINUE 
C  FINISHED. 

RETURN 

END 


DISTRIBUTION 


ADMINISTRATOR 

DEFENSE  DOCWENTATION  CENTER 
ATM  DDC-TCA  (12  COPIES) 

CAMERON  STATION,  BUILDING  S 
ALEXANDRIA,  VA  22314 

COMIANDER 

OS  ARMY  RSCH  t  STD  GP  (SUR) 

ATTN  LTC  JAMBS  N.  KENNEDY,  JR. 

CHIEF,  PHYSICS  G  NATH  BRANCH 
FPO  NW  YORK  09510 

COMMANDER 

OS  ARMY  ARMAMENT  MATERIEL 
READINESS  CCM1AMD 

ATTN  ORSAR-LIP-L,  TECHNICAL  LIBRARY 
ROCK  ISLAND,  IL  61299 

COMMANDER 

OS  ARMY  MISSILE  A  MUNITIONS 
CENTER  «  SCHOOL 
ATW  AT8K-CTD-F 
REDSTONE  ARSENAL,  AL  35809 

DIRECTOR 

OS  AMY  MATERIEL  SYSTEMS 
ANALYSIS  ACTIVITY 
ATM  DRXSY-MP 
ATTN  DRXSY-PO 

ABERDEEN  PROVING  GROUND,  MD  21005 
DIRECTOR 

US  AMY  BALLISTIC  RESEARCH  LABORATORY 
ATTN  DRDAR-T8B-8  (ST INFO) 

ATTN  DRXBR-AM,  W.  VAKAWTWERP 
ATTN  DR8TE-EL 
ATTN  CRQAR-BLB 

ABERDEEN  PROVING  GROUND,  W>  21005 

O. S.  ARMY  ELECTRONICS  TECHNOLOGY 
AND  DEVICES  IABORATORY 

ATTN  DELTT-OO 

FORT  MOMOVTH,  NJ  07703 

TEXAS  INSTRUMENTS,  INC. 

P. O.  BOX  226015 
ATTN  FRANK  POBLENZ, 

DALLAS,  TX  75266 

TBLEDYNE  BROUN  ENGINEERING 
CUMMINGS  RESEARCH  RARE 
ATTN  DR.  MELVIN  U  PRINCE,  MS-44 
HUNTSVILLE,  AL  35007 

ENGINEERING  SOCIETIES  LIEEARY 
345  EAST  47TB  STREET 
ATTN  ACQUISITIONS  DEPARTMENT 
MEN  YORK,  NY  10017 

DIRECTOR 

ARMED  POMCCS  RADIOBIOLOGY  RESEARCH 
INSTITUTE 

DEFENSE  NUCLEAR  AGENCY 
NATIONAL  NAVAL  MEDICAL  CENTER 
ATTN  RESEARCH  PROGRAM  COORDINATING 
OFFICER 

8ETHBSDA,  »  20014 


ASSISTANT  TO  THE  SECRETARY  OT  DBF  MSB 
ATOMIC  ENERGY 
ATTN  EXECUTIVE  ASSISTANT 
WASHINGTON,  DC  20301 

DIRECTOR 

DEFENSE  ADVANCED  RSCH  PROJ  AGENCY 
ARCHITECT  BUILDING 
ATTN  TIO 

1400  WILSON  BLVD. 

ARLINGTON,  VA  22209 

DIRECTOR 

DEFENSE  CIVIL  PREPAREDNESS  AGENCY 
ASSISTANT  DIRECTOR  FOR  RESEARCH 
ATTN  ADMIN  OFFICER 
ATTN  RE  (BO) 

ATTN  VO  (SB) 

WASHINGTON,  DC  20301 

DEFENSE  COMMUNICATIONS  ENGINEERING 
CENTER 

ATTN  CODE  R720,  C.  STANSBERRY 
ATTN  CODE  R123,  TECH  LIB 
ATTN  CODE  R400 
1860  W IE HIE  AVENUE 
RESTON,  VA  22090 

DIRECTOR 

DEFENSE  COMIUN IC ATION8  AGENCY 
ATTN  CCTC  C312 
ATTN  CODE  C313 
WASHINGTON,  DC  20305 

DIRECTOR 

DEFENSE  INTELLIGENCE  AGENCY 
ATW  RDS-3A 

ATTN  RD6-3A4,  POMPONIO  PLAZA 
WASHINGTON,  DC  20301 

DIRECTOR 

DEFENSE  NUCLEAR  AGENCY 

ATTN  RATH 

ATTN  DOST 

ATTN  RAEV 

ATTN  TITL 

ATTN  STVL 

ATTN  VLIS 

WASHINGTON,  X  20305 

COMIANDER 
FIELD  COMMAND 
DSFEN8E  NUCLEAR  AGENCY 
ATTN  FCPR 

ATTN  PCS  PM,  J.  SMITH 
ATTN  FCIMC 

KZRTLAND  AFB,  M  07115 
DIRECTOR 

IWrERSERVXCB  NUCLEAR  WEAPONS  SCHOOL 
ATTN  TTV 

KIRTLAND  AFE,  M  07115 

JOINT  CHIEFS  OF  STAFF 
ATTN  J-3 

WASHINGTON,  X  20301 


DIRECTOR 

JOINT  STRATEGIC  TARGET  PLANNING 
STAFF,  JCS 
OFFUTT  AFB 
ATTN  JSAS 
ATTN  JPST 

ATTN  NRI-STINPO  LIBRARY 
OMAHA,  m  60113 

CHIEF 

LIVERMORE  DIVISION 

FIELD  COMIAMD  DMA 

DEPARTMENT  OF  DEFENSE 

LAWRENCE  LIVEMORE  LABORATORY 

P.O.  BOX  000 

ATTN  FCPRL 

LIVEMORE,  CA  94550 

NATIONAL  COMIUNICATIONB  8YBTM 
OFFICE  OF  THE  MANAGER 
DEPARTMENT  OF  DEFENSE 
ATTN  HCS-TS,  CHARLES  D.  B0D60N 
WASHINGTON,  X  20305 

DIRECTOR 

NATIONAL  SECURITY  AGENCY 
DEPARTMENT  OF  DEFENSE 
ATTN  R-52,  O.  VAN  GUNTEN 
ATTN  S232,  D.  VINCENT 
FT.  MADE,  MD  20755 

UNDER  SBCY  OF  DBF  FOR  RSCH  6  ENGRG 
DEPARTMENT  OF  DEFENSE 
ATTN  G.  BARSE 
ATTN  StSS  (OS) 

WASHINGTON  DC  20301 

COMIANDER 
END  SYSTEM  OQN4AND 
DEPARTMENT  OF  THE  ARMY 
P.O.  BOX  1500 
ATTN  MDSC-AOLIB 
HUNTSVILLE,  AL  35807 

COMIANDER 

ERAXOM  TECHNICAL  SUPPORT  ACTIVITY 
DEPARTMENT  OF  THE  AMY 
ATTN  DELCS-K,  A  COHEN 
ATTN  DBLET-IR,  B.  HUNTER 
FORT  NORM  DOTH,  NJ  07703 

COMIANDER 

08  AMY  AMOR  CENTER 
ATTN  TECHNICAL  LIBRARY 
FORT  MOB,  KY  40121 

COMIANDER 

US  AMY  C0MM*CLBC  ENGRG  HMTAL 
AGENCY 

ATTN  COC-PRBO-S 
ATTN  00C-CBD-SE8 
FT  HUACBUCA,  AS  05613 

COMIAWDER 

00  AIVVY  COMMON  ZCATZOM8  OOMUUtD 
COMSAT  DEVELOPMENT  DIVISION 
ATTN  AT0I-CD-MO 
FT.  HUACNUCA,  AS  05613 


45 


DISTRIBUTION  (Cont'd) 


cany 

OB  AWT  CCMMDM 1CATIONS  STS  AGENCY 
ATTW  CCM-RD-T  COt-AD-SV 

fort  Monmouth,  mj  07703 
projnct  officer 

OS  AWT  COMMICATIOMS  ASS  A 
OSV  COMIAMD 
ATTN  mcm-ATC 
ATTN  DMCM-T08-BS1 
FORT  HOMOUTH,  NJ  07703 

DIVISION  ENGINEER 

IS  AWT  IHSIHtSR  D1V  HUNTS  VILLI 

P.O.  BOT  1600,  WIST  STATION 

ATTN  HNDED-SR 

ATTN  A.  T.  BOLT 

HUNTSVILLE,  AL  3S807 

OS  AWT  INTEL  THREAT  ANALYSIS 
OSTACIMERT 
ROOM  2201,  BLOG  A 
ARLINGTON  BALL  STATION 
ATTN  W  2200,  BLDG  A 
ARLINGTON,  VA  22212 

COMUNOKR 

IB  AWT  INTELLIGENCE  A  SIC  CKO 
ARLINGTON  HALL  STATION 
4000  ARLINGTON  BLVD 
ATTN  TECHNICAL  LIBRARY 
ATTN  TECH  INTO  FAC 
ARLINGTON,  TO  22212 

CaMUNDSR 

OS  ARMT  MISSILE  RESEARCH 

a  development  cam and 

ATTN  mcm-PB-EA,  WALLACE  O.  WAGNER 

ATTN  DNCW-FE-EG,  WILLIAM  S.  JOHNSON 

ATn  DREMI-TBD 

ATM  DRDMI-EAA 

REDSTONE  ARSENAL,  AL  3SS09 

OB  AWT  NUCLEAR  A  CHEMICAL  AGENCY 

7500  RAC  ALICE  ROAD 

BUILDING  2073 

ATTN  COL  A.  LOWRY 

AT«  DR.  J.  BERBERET 

SFRINGFICLD,  TO  22150 

CCSMAMDn 

OB  AWY  TEST  AMD  EVALUATION  COMMAND 
ATM  DRS1C-PA 

ABERDEEN  PROVING  GROUND,  MJ  21005 

COMMANDER 

OB  AWY  TRAINING  AMD 
OOCTRIMR  COMMAND 
ATW  ATORI-OP-8M 
FORT  MONROE,  TO  23651 

COMMANDER 

UNITE  SANDS  MISSILE  RANGE 
DEPARTMENT  OP  THE  AWY 
ATW  STEWO-W-AH,  J.  OEtMA 
NHITE  SANDS  MISSILE  RANGE,  W  MO 02 


OFF ZCER-IM -CHARGE 
CIVIL  ENGINEERING  LABORATORY 
WTOL  CONSTRUCTION  BATTALION  CENTER 
ATm  CODE  LOSA  (LIBRARY) 

ATTO  CODE  LOSA 

PORT  HUEMEKE,  CA  93041 


COMUNDER 

NAVAL  AIR  SYSTCMS  C0M1AMD 
ATn  AIR-350F 
WASHINGTON,  DC  21360 

COWANDER 

NAVAL  ELECTRCMIC  SYSTEMS  CCNIAMD 
at™  tm  117-215 
WASHINGTON,  DC  20360 

COMMANDER 

NAVAL  OCEAN  SYSTEMS  CENTER 
ATTO  COCK  015,  C.  FLETCHER 
ATW  RESEARCH  LIBRARY 
ATM  CODE  7240,  S.  N.  LICHTNAN 
SAN  DIEGO,  CA  92152 

CONUMDING  OFFICER 
NAVAL  ORDNANCE  STATION 
ATTN  STANDARDIZATION  DIV 
INDIAN  HEAD,  MD  20640 

SUPERINTENDENT  (CODE  1424) 

NAVAL  POSTGRADUATE  SCHOOL 
AT»  CODE  1424 
MONTEREY,  CA  93940 

DIRECTOR 

NAVAL  RESEARCH  LABORATORY 
ATTN  CODE  4104,  WAHUAL  L.  BRANCA  TO 
ATW  CODE  2627,  DORIS  R.  FOLEN 
ATM  CODE  6623,  RICHARD  L.  STATLER 
ATM  CODE  6624 
MASHDVOIOR,  DC  20375 

CORIANDER 

NAVAL  SNIP  ENGINEERING  CENTER 
DEPARTMENT  OP  THE  NAVY 
ATW  COBE  6174D2,  EDWARD  p,  DOTTY 
WASHINGTON,  DC  20362 

CORIANDER 

NAVAL  SURFACE  WEAPONS  CENTER 
ATM  CODE  F32,  EDWIN  R.  RATHBURM 
ATW  L.  LIBELLO,  CODE  WR43 
ATTN  CODE  NAS1RH,  W  130-100 
WHITE  OAK,  SILVER  SPRING,  MJ  20910 

CCmANDER 

WVAL  SURFACE  WEAPONS  CUTTER 
0AKLGRW  LABORAROTY 
ATW  CODE  DP-56 
DAELGREM,  TO  2244S 


CORIANDER 

NAVAL  WEAPONS  CENTER 

ATW  CODE  533,  TWCH  LIB 
CHINA  LAW,  CA  93555 


COMMANDING  OFFICER 

NAVAL  MAYORS  EVALUATION  FACILITY 

KIRTLAND  AIR  FORCE  EASE 

ATTO  CODE  AT -6 

ALBUQUERQUE,  W  87117 

OPTICS  OF  RATAL  RESRAACH 
ATTO  COBB  427 
ARLINGTON,  TO  22217 

DIRECTOR 

STRATEGIC  8TSTW8  PROJECT  OFFICE 
NAVY  DEPARTMENT 

attw  wp-2701,  John  w.  pitsenberger 
at™  NSP-2342,  RICHARD  L.  COLEMAN 
ATTO  NS  P-43,  TECH  LIB 
ATM  lisp-27334 
ATTO  NS  P-230,  D.  GOUJ 
WASHINGTON,  DC  20376 

C0HHAMDER 

AERONAUTICAL  SYSTEMS  DIVISION,  AP8C 
ATM  ASD-YR-EX 
ATTO  EHFTV 

WJUGHT -PATTERSON  APB,  OH  45333 

AIR  FORCE  TECHNICAL  APPLICATIONS 
CENTER 

ATM  TPS,  N.  SCHNEIDER 
PATRICK  APR,  PL  32925 

AF  WEAPONS  LABORATORY,  AFSC 
ATTN  NTH 
ATTO  NT 

Am  EL,  CARL  E.  BAM, 

ATTO  EL  XT 
ATTO  SUL 
ATW  CA 

ATTO  ELA,  J.  P.  CASTILLO 
ATTO  ELF 

ATTO  ELT,  W.  PAGE 
ATTO  MJB 

KIRTLAND  APB,  W  87117 
DIRSCTOR 

AIR  UNIVERSITY  LIBRARY 
department  of  the  air  ponce 
ATTO  AUL-LSE-7Q-2S0 
MAXWELL  AFB,  AL  36112 

WADQUANtBRS 

K1ECTROIIIC  SISTERS  DEVI 8 IOM/Y8EA 
WPARWENI  OP  TW  An  FORCE 
ATTO  YSEA 

HARBCQM  APR,  NA  01731 
COMARDER 

FOREIGN  TECHNO  I/3GY  DIVISION,  AFSC 
ATTO  NICD  LIBRARY 
ATW  STOP,  S.  U  BALLARD 
NR1GNT-PRTTERSQN  APB,  OR  45433 

COWtARDER 

OGDEN  ALC/MNEDOE 

DEPARMENI  or  TW  An  PONCE 

ATTO  OO-ALC/WHTH,  P.  W.  BERTNEL 

ATW  NNEDO,  LEO  KIDMAN 

ATM  NAJ  R.  BMCKBUW 

RILL  ATI,  UT  444 06 


46 


DISTRIBUTION  (Cont'd) 


CGMUNOSK 

rone  air  development  c  wm  at  sc 

ATTM  TSLD 

(Mrrriss  kn,  mr  13441 

CQWUkMMR 

RACRAMEVTO  AIR  L0GX8TXC8  CUTTER 
DEPARTMENT  OR  TB  AIR  FONCE 
ATTN  NNCRS,  R.  A.  PEUIA8TRO 
ATTN  NMIRA,  J.  V.  DNMES 

attn  msR m,  r.  r.  spear 

NOCLRLLMI  AFB,  CA  95652 
EAMSO/IH 

AIR  FORCE  SYSTEMS  COMMAND 
PONT  OFFICE  RON  92960 
MORLDVAY  POSTAL  CENTER 
( INTELLIGENCE) 

ATTN  DID 

LOS  ANGELES,  CA  90009 
8AM80/MN 

AIR  FORCE  SYSTEMS  C0M1AND 
(MWOTEMAN) 

ATTN  MNNH,  MAJ  N.  SARAH 
ATTN  MNNH,  CAPT  R.  I.  LAWRENCE 
NORTON  APB,  CA  92409 

SAMSO/YA 

AIR  FORCE  SYSTEMS  COMMAND 
POST  OFFICE  B0K  92960 
WORLDNAY  POSTAL  CENTER 
ATTN  YAPC 

LOS  ANGELES,  CA  90009 

STRATEGIC  AIR  COMAND/XFFS 
ATTN  NRI-STINFO  LIBRARY 
ATTN  DEL 

ATTN  GARNET  E.  MATZKE 

ATTN  XPF8,  MAJ  BRIAN  G.  STEPHEN 

OTFUTT  APB,  MB  68113 

DEPARTMENT  0T  ENERGY 
ALBUQUERQUE  OPERATIONS  OPFICE 
P.O.  BOX  5400 

ATTN  DOC  CON  FOR  TECH  LIBRARY 
ATTN  OPERATIONAL  SAFETY  DIV 
ALBUQUERQUE,  W  07t15 

UNIVERSITY  OF  CALIFORNIA 
LAWRENCE  LIVERMORE  LABORATORY 
P.O.  BOX  SOB 

ATTN  DOC  CON  FOR  TECHNICAL 
INFO  WAT  ION  DEPT 

ATTN  DOC  CON  FOR  L-06,  T.  DONICH 
ATTN  DOC  CON  FOR  L-545,  D.  MEEKER 
ATTN  DOC  CON  FOR  L-156,  E.  MILLER 
ATTN  DOC  CON  FOR  L-10,  H.  KRUGER 
ATTN  DOC  CON  FOR  H.  8.  CABAYAN 
LIVERMORE,  CA  94550 

LOS  AMMO*  SCIENTIFIC  LABORATORY 
P.O.  SOK  1663 

ATTN  DOC  CON  FOR  BRUCE  N.  NOEL 
ATTN  DOC  CON  FOR  CLARENCE  BENTON 
LOS  AIAM08,  m  87545 


k — i  ■  -  — »  -  -  • — ■-  ■  iMk am  mti’r 


SAMDIA  LABORATORIES 
P.O.  BOX  5800 

ATTN  DON  CON  FOR  C.  N.  VXTTXTOB 
ATTN  DON  CON  POR  R.  L.  PARKER 
ATTN  DOC  CON  POR  EZJtER  P.  HARTMAN 
ALBUQUERQUE,  W  87115 

CENTRAL  INTELLIGENCE  AGENCY 
ATTN  RD/SI,  Ml  5G40,  HQ  BLDG 
FOR  06I/MSD/MNB 
WASHINGTON,  DC  20  505 

ADMINISTRATOR 

DEFEN8E  ELECTRIC  POWER  ADMIN 
DEPARTMENT  OF  THE  INTERIOR 
INTERIOR  SOUTH  BLDG,  312 
ATTN  L.  O'NEILL 
WASHINGTON,  DC  20240 

DEPARTMENT  GF  TRANSPORTATION 
FEDERAL  AVIATION  ADMINISTRATION 
HEADQUARTERS  SBC  DIV,  ASX-300 
800  INDEPENDENCE  AVENUE,  SW 
ATTN  SC  DIV  ASE-300 
WASHINGTON,  DC  20591 

AEROSPACE  CORPORATION 
P.O.  BOX  92957 
ATTN  C.  B.  P1ARLSTON 
ATTN  IRVING  M.  GARFUNKBL 
ATTN  JULIAN  RBINHBIMXR 
ATTN  LIBRARY 
ATTN  CHARLES  GREENHOW 
LOS  ANGELES,  CA  90009 

AGBABIAN  ASSOCIATES 
250  NORTH  MASH  STREET 
ATTN  LIBRARY 
EL  8BGUNDO,  CA  90245 

AVCO  RESEARCH  *  SYSTWS  GROUP 
201  LOWELL  STREET 
ATTN  N.  LEPSEVICH 
WIWINOPON,  HA  01887 

BATTELXB  MEMORIAL  INSTITUTE 
505  KING  AVENUE 
ATTN  ROBERT  H.  BALZEK 
ATTN  EUGENE  R.  IEACH 
COLUMBUS,  OH  43201 

BDM  CORPORATION 
7915  JONES  BRANCH  DRIVE 
ATTN  CORPORATE  LIBRARY 
MCLEAN,  VA  22101 

BOM  CORPORATION 
P.O.  BOX  9274 
ALBUQUERQUE  INTERNATIONAL 
ATTN  LIB 

ALBUQUERQUE,  NN  87119 

BENDIX  CORPORATION,  THE 
RESEARCH  LABORATORIES  DIVISION 
BENDIX  CENTER 
ATTN  MAX  FRANK 
SOUTHFIELD,  HI  48075 


BENDIX  CORPORATION 
NAVIGATION  AND  CONTROL  GROUP 
ATTN  DBPT  6401 
TBTBRBORO,  NJ  07608 

BONING  COMPANY 
P.O.  BOX  3707 
ATTN  HOWARD  N.  IfXCKLXIM 
ATTM  D.  E.  ISBELL 
ATTN  DAVID  KWLB 
ATTM  B.  C.  HANRABAM 
ATTN  KENT  TECH  LIB 
8 BATTLE,  HA  98124 

B008-ALLEH  AMD  HAMILTON,  INC. 

106  APPIE  STREET 
ATTM.  JU.  J.  CHRX8NER 
ATM  TECH  LIB 
TXMTOM  FALLS,  NJ  07724 

BURROUGHS  CORPORATION 
FEDERAL  AMD  SPECIAL  SYSTSM8  GROUP 
CENTRAL  AVE  AMD  ROOT*  252 
P.O.  BOX  517 

ATTN  ANGELO  J.  NAURIELLO 
PAOLX,  PA  19301 

CALS  PAN  CORPORATION 
P.O.  BOX  400 
ATTM  TECH  LIBRARY 
BUFFALO,  NY  14225 

CHARLES  STARK  DRAPER  LABORATORY  INC. 

555  TECHNOLOGY  SQUARE 

ATTN  KENNETH  FERTIG 

ATTM  TIC  MS  74 

CAMBRIDGE,  MA  02139 

CINCINNATI  ELECTRONICS  CORPORATION 
2630  GLENDALE-NXLFORD  ROAD 
ATTN  LOIS  HAMtOHD 
CINCINNATI,  OH  45241 

COMPUTER  SCIENCES  CORPORATION 
6565  ARLINGTON  BLVD 
ATTN  RAMONA  BRIGGS 
FALLS  CHURCH,  VA  22046 

COMPUTER  SCIENCES  CORPORATION 
1400  SAN  MATEO  ELVD,  SB 
ATTN  RICHARD  B.  DXCKHAOT 
ATTN  ALVIN  SCHIFF 
ALBUQUERQUE,  Ml  87108 

CONTROL  DATA  CORPORATION 
P.O.  SOX  O 
ATTN  JACK  MEEHAN 
MINNEAPOLIS,  NN  55440 

CUrtER-HANMER,  INC. 

AIL  DIVISION 
CQNAC  ROAD 
ATTN  EONAW  KARPIN 
DEER  PARK,  MY  11729 

DIKBNOOO  INDUSTRIES,  INC 
1009  BRANOBURY  DRIVE,  SB 
ATTN  TECH  LIB 
ATTN  L.  WAYNE  DAVIS 
ALBUQUERQUE,  IM  87108 


47 


.Auatjtw.  affMiMfiir*  .y  nil  irii# 


DISTRIBUTION  (Cont'd) 


dxkewooo  industries,  me. 

1100  GXJtNlXtt  MW 01 

ATTN  K.  US 

LOB  ANGELES,  CA  00024 

K-SYSTM8,  INC 
GREENVILLE  DIVISION 
9.0.  BOX  1056 
ATTN  JOLBTA  NOONS 
GREENVILLE,  TX  75401 

BFPKTS  TECHNOLOGY.  INC. 

5383  HOLLISTSS  AVSMUl 

ATTN  S.  CLOU 

SANTA  BARBARA,  CA  93111 

EGAG  WASHINGTON  ANALYTICAL 
SERVICES  CENTER,  INC. 

P.O.  BOX  10218 
ATTN  C.  GILES 
ALBUQUERQUE,  m  87114 

EXXON  NUCLEAR  COMPANY,  INC. 
RESEARCH  AND  TECHNOLOGY  CENTER 
2955  GEORGE  WASHINGTON  WAY 
ATTN  DR,  A.  W,  TRXVtLPXBCS 
RICHLAND,  WA  99352 

FAIRCHILD  CAMERA  AMD  INSTRUMENT 
464  ELLIS  STREET 
ATTN  SBC  CON  FOR  DAVID  R.  MYERS 
MOUNTAIN  VIEW,  CA  94040 

FORD  AEROSPACE  6  COMMUNICATIONS 
3939  FABIAN  WAY 
ATTN  TECHNICAL  LIBRARY 
PALO  AUTO,  CA  94303 

FORD  AEROSPACE  4  COMMUNICATIONS 
OPERATIONS 

FORD  6  JAMBOREE  ROADS 
ATTN  KEN  C.  ATTINGER 
ATTN  B.  R.  PONCELET,  JR. 

NWPORT  BEACH,  CA  92663 

FRANKLIN  INSTITUTE,  THE 
20TH  STREET  AND  PARKWAY 
ATTN  RAMIE  H.  THOMPSON 
PHILADELPHIA,  PA  19103 

GENERAL  DYNAMIC 8  CORF 
ELECTRONICS  DIVISION 
P.O.  BOX  SI 125 
ATTN  RECK  LIB 
SAN  DIEGO,  CA  92138 

GENERAL  DYNAMICS  CORPORATION 

INTER-DIVISION  RESEARCH  LIBRARY 

KEARNY  MESA 

P.O.  BOX  00847 

ATTN  RESEARCH  LIBRARY 

SAN  DIEGO,  CA  96123 

GENERAL  ELECTRIC  CO.-TWFO 
CENTER  FOR  ADVANCED  STUDIES 
816  8TATE  STREET  <PO  DRAWER  QQ) 
ATTN  QA81AC 

ATTN  HOYDEN  A.  RUTHERFORD 
ATTN  WILLIAM  NCNAMERA 
SANTA  BARBARA,  CA  93102 


GENERAL  ELECTRIC  COMPANY 
AEROSPACE  ELECTRONICS  SY8TSIS 
FRENCH  HOAD 

ATTN  CHARLES  M.  HEWXSON 
UTICA,  NY  13503 

GENERAL  ELECTRIC  COMPANY 
P.O.  BOX  $000 
ATTN  TECH  LIB 
BINGHAMTON,  NY  13902 

GENERAL  ELECTRIC  CO. -TEMPO 
ALEXANDRIA  OFFICE 
HUNTINGTON  BUILDING,  SUITE  300 
2560  HUNTINGTON  AVENUE 
ATTN  DAS I AC 
ALEXANDRIA,  VA  22303 

GENERAL  RESEARCH  CORPORATION 

SANTA  BARBARA 

P.O.  BOX  6770 

ATTN  TECH  INFO  OFFICE 

SANTA  BARBARA,  CA  93111 

GEORGIA  INSTITUTE  OF  TECHNOLOGY 
GEORGIA  TECH  RESEARCH  INSTITUTE 
ATTN  R.  CURRY 
ATLANTA,  GA  30332 

CORP 

(SORGIA  INSTITUTE  OF  TECHNOLOGY 
OFFICE  OF  CONTRACT  ADMINISTRATION 
ATTN  RES  4  SBC  COORD  FOR  HUGH  DENNY 
ATLANTA,  GA  30332 

CORP 

GRtPMAN  AEROSPACE  CORPORATION 
SOUTH  OYSTER  BAY  ROAD 
ATTN  L-01  35 
BBTHPAGE,  NY  11714 

GTE  SYLVANIA  INC. 

ELECTRONICS  SYSTEMS  GRP-EA STERN  DIV 
77  A  STREET 

ATTN  CHARLES  A.  THORNHILL,  LIBRARIAN 
ATTN  LBOHDARD  L.  BLAIS DELL 
NEEDHAM,  HA  02194 

GTE  SYLVANIA,  INC. 

189  8  STREET 

ATTN  CHARLES  H,  RAMSBOTTGM 

ATTN  DAVID  D.  FLOOD 

ATTN  BfXL  P,  MOTCHOK 

ATTN  H  4  V  GROUP,  MARIO  A.  NUREPORA 

ATTN  J.  WALDRON 

NEEDHAM  HEIGHTS,  MA  02194 

KARRIS  CORPORATION 

KARRIS  8D4ICONDUCTOR  DIVISION 

P.O.  BOX  883 

ATTN  V  PRES  4  MGR  PROMS  DIV 
MELBOURNE,  PL  32901 

HAZELTINE  CORPORATION 
PULA8KI  ROAD 

ATTN  TECH  INFO  CTR,  M.  WAITE 
GREEN LAWN,  NY  11740 


HONEYWELL  INCORPORATED 
AVIONICS  DIVISION 
2600  RIDGEWAY  PARKWAY 
ATTN  84  RC  LIB 
ATTN  RONALD  R.  JOHNSON 
MINNEAPOLIS,  MN  55413 

HONEYWELL  INCORPORATED 

AVIONICS  DIV 18  ION 

13350  U.8.  HIGHWAY  19  WORTH 

ATOI  M.S  725-5,  STACEY  H.  GRAFF 

ATTN  W.  K  STEWART 

ST.  PETERSBURG,  FL  33733 

HUGHES  AIRCRAFT  COMPANY 
CBWTINBLA  AND  TBALE 
ATTN  Ctom  B.  SINGLETARY 
ATTN  CTDC  6/El 10 
ATTN  KENNETH  R.  WALKER 
CULVER  CITY,  CA  90230 

I  IT  RESEARCH  INSTITUTE 
ELECTROMAG  C0MPATA8ILITY  ANAL  CTR 
NORTH  SEVERN 
ATTN  RCOAT 
ANNAPOLIS,  MD  21402 

I IT  RESEARCH  INSTITUTE 
10  WEST  35TH  STREET 
ATTN  IRVING  N.  MIN DEL 
ATTN  JACK  E.  BRIDGES 
CHICAGO,  XL  60616 

INSTITUTE  FOR  DEFENSE  ANALYSES 
400  ARMY-NAVY  DRIVE 
ATTN  TECH  INFO  8ERVICES 
ARLINGTON,  VA  22202 

INTL  TEL  4  TELEGRAPH  CORPORATION 
500  WASHINGTON  AVENUE 
ATTN  TECHNICAL  LIBRARY 
ATTN  ALEXANDER  T.  RICHARDSON 
NUT LEY,  NJ  07110 

IRT  CORPORATION 
P.O.  BOX  81087 
ATTN  C.  B.  WILLIAMS 
ATTN  DENNIS  SWIFT 
SAN  DIEGO,  CA  92138 

JAYCOR 

SANTA  BARBARA  FACILITY 
P.O.  BOX  2008 
ATTN  W.  A.  RADA SKY 
SANTA  BARBARA,  CA  93120 

JAYCOR 

1401  CAM INO  DEL  MAR 
ATTN  ERIC  P.  WENAAS 
ATTN  RALPH  H.  STAHL 
DSL  MAR,  CA  92014 

JAYCOR 

205  S  WHITING  STREET,  SUITS  500 
ATIM  LIB 

ALEXANDRIA,  VA  22304 


DISTRIBUTION  (Cont'd) 


■UMAX  SCIENCES  CORPORATION 
1500  GARDE*  OF  TR  0006  ROAD 
ATTN  AIR  CRT  P.  BRIDGES 

Atm  v-  roam  rich 

ATTN  HALTS*  E.  HARE 
ATTH  PRARK  H.  SHELTON 
ATT*  JERRY  I.  DWELL 
ATT*  PHIL  TRACY 
ATT*  HERMER  STARK 
COLORADO  SPRINGS,  CO  SOSO 7 

LITTON  SYSTEMS,  IMC. 

DATA  SYSTEMS  DIVISION 
9000  WOODLEY  AVENUE 
ATTN  WC  GP 
ATTN  M848-61 
VAN  NUTS,  CA  91409 

LITTON  SYSTEMS,  INC. 

AMECON  DIVISION 
5115  CALVERT  ROAD 
ATTN  J.  SKAGGS 
COLLEGE  PARK,  HD  20740 

r^wrimn  MISSLES  AMD  SPACE 
COMPANY,  INC 
P.O.  BOX  504 
ATTN  L.  ROSSI 
ATTN  SAMUEL  I.  TAIMITPY 
ATTN  H.  E.  THAYN 
ATTN  GEORGE  P.  HEATH 
ATTN  BENJAMIN  T.  KIMURA 
SUNNYVALE,  CA  94086 

LOCKHEED  MISSILES  AND  SPACE 
COMPANY,  INC. 

3251  HANOVER  STREET 
ATTN  TECH  INTO  CTR  D/COLL 
PALO  ALTO.  CA  94304 

M.I.T.  LINCOIN  LABORATORY 
P.O.  BOX  73 
ATTN  LEONA  LOUGHLIN 
LEXINGTON,  MA  02173 

MARTIN  MARIETTA  CORPORATION 
ORLANDO  DIVISION 

P.O.  BOH  5837 

ATTN  MONA  C.  GRIFFITH 

ORLANDO,  PL  32805 

MCDONNEL  DOUGLAS  CORPORATION 
POST  OFFICE  BOX  516 
ATTN  TOM  ENOCH 
ST.  LOUIS,  MO  63166 

MCDONNELL  DOUGLAS  CORPORATION 
5301  BOLSA  AVENUE 
ATTN  STANLEY  SCHNEIDER 
ATTN  TECH  LIBRARY  SERVICES 
HUNTINGTON  BEACH,  CA  92647 

MISSION  RESEARCH  CORPORATION 

P.O.  DRAWER  719 

ATTN  W  GROUP 

ATTN  WILLIAM  C.  WRT 

ATTN  C.  LONOIIRE 

SANTA  BARBARA,  CA  93102 


MISSION  RESEARCH  CORPORATION 
m  SYSTB1  applications  division 
1400  SAN  MATEO  BLVD,  SE,  SUITE  A 
ATTN  DAVID  E.  MERIWETHER 
ATTN  L.  N.  MCCORMICK 
ALBUQUERQUE,  181  87108 

MISSION  RESEARCH  CORPORATION 
-SAN  DIEGO 
P.O.  BOX  1209 
ATTN  V.  A.  J.  VAN  LINT 
IA  JOLLA ,  CA  92038 

MITRE  CORPORATION,  THE 
P.O.  BOX  208 
ATTN  M.  P.  FITZGERALD 
BEDFORD,  MA  01730 

NORDEN  SYSTEMS,  INC. 

HELEN  STREET 

ATTN  TECHNICAL  LIBRARY 

NORWALK,  CT  06856 

NORTHROP  RESEARCH  TECHNOLOGY  CENTER 
ONE  RESEARCH  PARE 
ATTN  LIBRARY 

PALOS  VERDES  PENN,  CA  90274 

NORTHROP  CONTORTION 
ELECTRONIC  DIVISION 
2301  NEST  120TH  STREET 
ATTN  LEV  SMITH 
ATTN  RAD  EFFECTS  GRP 
HAWTHORNS,  CA  90250 

PHYSICS  INTERNATIONAL  COMPANY 
2700  MERCED  STREET 
ATTN  DOC  CON 
SAN  UAHDRO,  CA  94577 

R.t.O  ASSOCIATED 

P.O.  BOX  9695 

ATTN  S.  CLAY  ROGERS 

ATTN  RICHARD  R.  SCHAEFER 

ATTN  DOC  CON 

ATTN  M.  GROVER 

ATTN  C.  MACDONALD 

MARINA  DEL  HEX,  CA  90291 

RAD  ASSOCIATES 
1401  WILSON  BLVD 
SUITE  500 
ATTN  J.  BOMBARDS 
ARLINGTON,  VA  22209 

RAND  CORPORATION 
1700  NAM  STREET 
ATTN  LIE-0 
ATTN  N.  BO  LLP  HEY 
SANTA  MONICA,  CA  90406 

RAYTHEON  COMPANY 
HARWELL  MAD 
ATT*  OA JAMAH  W.  JOEI 
REDFOND,  MA  01730 

AAYTWOM  CWPAMY 
528  BOSTON  POST  SOAP 
ATT*  REMOLD  U  FLEECE* 

SUDBURY,  HA  0177* 


RCA  CORPORATION 

GOVERNMENT  SYSTBtS  DIVISION 

ASTRO  ELECTRONICS 

P.O.  BOX  800,  LOCUST  CORNER 

EAST  WINDSOR  TOWNSHIP 

PRINCETON,  NJ  08540 

RCA  CORPORATION 

DAVID  SARMOFF  RESEARCH  CENTER 

P.O.  BOX  432 

ATTN  SECURITY  DEPT,  L.  WHICH 
PRINCETON,  NJ  08540 

RCA  CORPORATION 
CAMDEN  COMPLEX 
PROMT  6  COOPER  STREETS 
ATTN  OLIVE  WHITEHEAD 
ATTN  R.  N.  TOSTROM 
CAMDEN,  NJ  08012 

ROCKWELL  INTERNATIONAL  CORPORATION 

P.O.  BOX  3105 

ATTN  N.  J.  RUDIE 

ATTN  J.  L.  MONROE 

ATTN  V.  J.  MICHEL 

ATTN  D/243-068,  031-CA31 

ANAHEIM,  CA  92803 

ROCKWELL  INTERNATIONAL  CORPORATION 
SPACE  DIVISION 

12214  SOUTH  LAKEWOOD  BOULEVARD 
ATTN  B.  E.  WHITE 
DOWNEY,  CA  90241 

ROCKWELL  INTERNATIONAL  CORPORATION 

815  LAPAIK  STREET 

ATTN  B— 1.  DIV  TIC  (BAGS) 

EL  SEGUNDO,  CA  90245 

ROCKWELL  INTERNATIONAL  CORPORATION 
P.O.  BOX  369 
ATTN  P.  A.  SHAH 
CLEARFIELD,  UP  84015 

SANDERS  ASSOCIATES,  INC. 

95  CANAL  STREET 

ATTN  1-6270,  R.  G.  DBS  PATH! .  SR  P  E 
NASHUA,  MH  03060 

SCIENCE  APPLICATIONS,  INC. 

P.O.  BOS  277 

ATTN  FREDERICK  N.  TE9CHK 

BERKELEY,  CA  94701 

SCIENCE  APPLICATIONS,  INC. 

P.O.  ROX  2351 
ATTN  R.  PARKINSON 
LA  JOLLA,  CA  9203E 

SCIENCE  APPLICATIONS,  INC. 
HUNTSVILLE  DIVISION 
2109  N.  CLINTON  AVENUE 
SUITE  700 
ATTN  NOEL  E.  BYRE 
HUNTSVILLE,  AL  35805 

ICIERCE  APPLICATIONS,  INC. 

MOO  VESTPAAX  DRIVE 
ATTN  WILLIAM  U  CHADS EY 
MCLEAN,  VA  22101 


49 


DISTRIBUTION  (Cont’d) 


singer  company 

ATTN:  SECURITY  MANAGER 
FOR  TECH  INFO  CTR 
1150  HC  BRIO*  AVENUE 
LITTLE  FALLS,  RJ  07424 

SPERRY  RAND  CORPORATION 
SPERRY  MICROWAVE  ELECTRONICS 
P.O.  BOX  4648 
ATTM  MARGARET  COPT 
CLEAWATXR,  FL  33518 

SPERRY  RAMD  CORPORATION 
SPERRY  DIVISION 
MARCOS  AVENUE 
ATTN  TECS  LIB 
GREAT  NECK,  NY  11020 

SPERRY  RAMD  CORPORATION 
SPERRY  FLIGHT  SYSTEMS 
P.O.  BCE  21111 
ATTN  D.  ANDREW  SCHW 
PHOENIX,  AZ  85036 

SPIRE  CORPORATION 
P.O.  BOX  D 
ATTN  JOHN  R.  UGLUH 
ATTN  ROGER  G.  LITTLE 
BEDFORD,  NA  01730 

SRI  INTERNATIONAL 
333  RAVENBMOOD  AVENUE 
ATTN  ARTHUR  LEE  WHITSON 
MENLO  PARE,  CA  94025 

SYSTEMS,  SCIENCE  AMD  SOFTWARE,  INC. 
F.O.  BOX  1620 
ATTN  ANDRM  R.  WILSON 
LA  JOLLA,  CA  92038 

TEXAS  INSTRUMENTS,  INC. 

P.O.  BOX  6015 
ATTN  TECH  LIB 
ATTN  DONALD  J.  MANUS 
DALLAS.  TX  75265 

TEN  DEFENSE  A  SPACE  SYS  GROUP 

ONE  SPACE  PANE 

ATTN  O.  E.  ADAMS 

ATYN  R.  L  PLEBUCH 

ATTN  L.  R.  MAGNOLIA 

ATTN  1L  R.  HOLLOWAY 

ATTN  W.  a ABOARD 

REDONDO  BEACH,  CA  90278 


TEXAS  TECH  UNIVERSITY 

P.O.  BOX  5404  NORTH  COLLEGE  STATION 

ATTN  TRAVIS  L.  SIMPSON 

LOBBOCX,  TX  79417 

UNITED  TECHNOLOGIES  CORP 
HAMILTON  STANDARD  DIVISION 
BRADLEY  INTERNATIONAL  AIRPORT 
ATTN  CHIEF  ELHC  DESIGN 
WINDSOR  LOCKS ,  CT  06069 

WESTINGNOUSE  ELECTRIC  CORPORATION 
ADVANCED  ENERGY  SYSTEMS  DIV 
P.O.  BOX  10B64 
ATTM  TECH  LIB 
PITTSBURGH,  PA  15236 

US  AMY  ELECTRONICS  RESEARCH 
A  DEVELOPMENT  COWARD 
ATTN  TECHNICAL  DIRECTOR,  DRDRL-CT 

HARRY  DIAMOND  LABORATORIES 

ATTN  CO/TD/TSO/DIVI8ION  DIRECTORS 

ATTN  RECORD  COPY,  B1200 

ATTN  HDL  LIBRARY  81100  (3  COPIES) 

ATTN  HDL  LIBRARY  (NOODBRIDGE)  81100 

ATTN  TECHNICAL  REPORTS  BRANCH,  81300 

ATTN  CHAX1BIAN,  EDITORIAL  COMIITTEE 

ATTM  CHIEF,  21000 

ATTN  CHIEF,  22000 

ATTN  CHIEF,  22100  (3  COPIES) 

ATTN  CHIEF,  22300 
ATTN  CHIEF,  22800 
ATTN  CHIEF,  22900 
ATTM  CHIEF,  13300 
ATTN  CHIEF,  21100  (3  COPIES) 

ATTM  CHXHF,  21200 

ATTN  CHIEF,  21300  (5  COPIES) 

ATTN  CHIEF,  21400  (2  COPIES) 

ATTN  CHIEF,  21500 

ATTN  BALXCKI,  F.  N.,  20240 

ATTN  NDIXNITZ,  F.  H.,  20240 

ATTN  TALLEHICO,  A.,  47400 

ATTN  BIXBY,  N.,  22900 

ATTN  WYATT,  N.  T.,  21300  (20  COPIES) 


50 


